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Abstract 

Mutations and oxidative modification in the protein Cu,Zn superoxide dismutase (SOD1) have been implicated in the death 
of motor neurons in amyotrophic lateral sclerosis (ALS), a presently incurable, invariably fatal neurodegenerative disease. 
Here we employ steered, all-atom molecular dynamics simulations in implicit solvent to investigate the significance of either 
mutations or post-translational modifications (PTMs) to SOD1 on metal affinity, dimer stability, and mechanical malleability. 
The work required to induce moderate structural deformations as a function of sequence index constitutes a "mechanical 
fingerprint" measuring structural rigidity in the native basin, from which we are able to unambiguously distinguish wild- 
type (WT) SOD1 from PTM variants, and measure the severity of a given PTM on structural integrity. The cumulative 
distribution of work values provided a way to cleanly discriminate between SOD1 variants. Disulfide reduction destabilizes 
dimer stability more than the removal of either metal, but not moreso than the removal of both metals. Intriguingly, 
we found that disulfide reduction mechanically stabilizes apo SOD1 monomer, underscoring the differences between native 
basin mechanical properties and equilibrium thermodynamic stabilities, and elucidating the presence of internal stress in the 
apo state. All PTMs and ALS-associated mutants studied showed an increased tendency to lose either Cu or Zn, and to 
monomerize- processes known to be critical in the progression of ALS. The valence of Cu strongly modulates its binding free 
energy. As well, several mutants were more susceptible to loss of metals and monomerization than the disulfide-reduced or 
apo forms of SOD1. Distance constraints are required to calculate free energies for metal binding and dimer separation, which 
are validated using thermodynamic cycles. When distance constraints are removed, the results agree with those obtained 
from direct application of the Jarzynski equality. 
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1 Introduction 



Copper, zinc superoxide dismutase (SOD1) is a homo- 
dimeric antioxidant enzyme of 32kDa present in all eukary- 
otes. Each monomer consists of an eight stranded greek key 
(3 barrel of 153 amino acids [1-3] and binds one Cu and one 
Zn ion, which significantly enhance native thermodynamic 
and mechanical stability [4,5]. 

Two large, functionally important loops determine the struc- 
ture and activity of the enzyme. Loop VII or the electro- 
static loop (residues 121-142) contains charged and polar 
residues that enhance enzymatic activity by inducing an 
electrostatic funnel towards the active site centered on the 
Cu ion [6], which catalyzes the conversion of superoxide 
(0^~) to less toxic species(either O2 or H2O2). Loop IV or 
the Zn-binding loop (residues 49-83) coordinates Zn with 
histidines 63, 71, 80 and Aspartic acid 83, which, along with 
a disulfide bond between Cysteines 57 and 146, enforce con- 
comitant tertiary structure in the protein [7]. Residues in 
the Zn binding loop form about 38% of the dimer surface 
contact area, so that disorder in the loop due to Zn expul- 
sion and/or disulfide reduction can facilitate monomeriza- 
tion of the homodimer [8-12]. 

The Cu on the other hand is coordinated by residues mainly 



in the (3 strands of the immunoglobulin-like core of the 
protein- histidines 46, 48, 120, along His 63 which bridges 
the Cu and Zn ions- so that protein structure only weakly 
couples to Cu binding [13,14], supporting a primarily enzy- 
matic role of Cu. 



Amino acid missense mutations at more than 150 positions 
in SOD1 have been found to cause amyotrophic lateral scle- 
rosis (ALS), an invariably fatal neurodegenerative disease 
characterized by loss of the motor neurons in the brain, 
brainstem and spinal cord (http://alsod.iop.kcl.ac.uk/). 
Such familial mutations constitute about 20% of the cases 
of ALS known to display autosomal dominance [15], con- 
ferring symptoms through a mechanism involving a toxic 
gain of function, in that SOD1 knockout mice do not de- 
velop the neurodegeneration indicative of the familial ALS 
(fALS)-like phenotype [16]. Gain of function symptoms have 
been variously attributed to generation of reactive oxygen 
and nitrogen species, cytoskeletal disruption, caspase acti- 
vation, mitochondrial dysfunction, proteosome disruption, 
and microglial activation [17-19]. SOD-mediated fALS is 
one of the most prominent identified causes of the disease 
despite constituting only approximately 2-5% of all known 
cases. The vast majority, about 90%, of ALS cases have no 
known underlying etiology and are termed sporadic (sALS). 
Macroscopically, sALS is clinically indistinguishable from 
fALS [20]. Lewy body-like inclusions in sALS have been 
found to be immunoreactive to SODl-specific antibod- 



ies [21], and misfolding-specific antibodies have also iden- 
tified misfolded SOD1 protein in sALS inclusions [22,23]. 
At the molecular level however, protease resistant cores 
in aggregates as identified by MALDI-TOF mass analysis 
were observed to differ between the WT sequence and the 
fALS-associated mutants G37R, G85R, and G93A [24], 
indicating that mutations can modulate structural poly- 
morphism in aggregate cores, resulting in distinct physico- 
chemical aggregate properties such as sequence-dependent 
solubility. The mutation-induced polymorphism of fibril 
cores has been recapitulated in coarse-grained simulations 
using discrete molecular dynamics with structure-based, 
Go-like potentials [25] Nevertheless, several studies have 
indicated oxidative modification of WT SOD1 can be toxic 
in ALS, linking SOD1 misfolding to a common pathogenic 
mechanism in fALS and sALS [26-30]. 

Studies involving misfolding-specific antibodies have shown 
that misfolded SOD1 (either G85R or G127X) can in- 
duce misfolding in natively-structured WT SOD1 by direct 
protein-protein interactions [31]. Holo, pseudo WT SOD1 
(C6A/C111S) has been observed to aggregate under phys- 
iologically relevant conditions, with characteristics similar 
to aggregates in fALS patients [32]. Further, cell to cell 
transfer mediated by macropinocytosis has been observed 
in the fALS mutant H46R [33]. These studies point to pri- 
onogenic mechanisms for the propagation of SOD1 misfold- 
ing and aggregation- a common theme in protein-misfolding 
diseases. 

SOD1 protein need not globally misfold or unfold in or- 
der to aggregate. Accumulated evidence for several pro- 
teins indicates that aggregation may be initiated from lo- 
cally (rather than globally) unfolded states [34]. These par- 
tially disordered states may be induced by external agents 
or may become accessible via thermal fluctuations or rare 
events. For the case of SOD1, a near native aggregation 
precursor was found for an obligate monomeric SOD1 vari- 
ant (C6A/C111A/F50E/G51E), wherein protective cap re- 
gions are locally unfolded around the native /3-stranded 
core [35]. As well, the crystal structures of the fALS mu- 
tants S134N and apo H46R have both been observed to 
adopt fibrillar structures with intercalating loops between 
largely native-like domains [36], and S134N has been ob- 
served to form oligomers in solution that are stabilized by 
locally unfolded elements from the native structure, involv- 
ing transient interactions between electrostatic loops from 
different dimers [37]. 

The above studies have motivated a previous computa- 
tional study, where we had focused on the native and near- 
native mechanical properties of WT SOD1, premature vari- 
ants lacking post-translational modifications, and both ALS- 
associated and rationally designed truncation mutants [5]. 
We felt that this approach was complimentary to experi- 
mental structural and thermodynamic assays. In that work, 
we examined relatively large changes in native structure as 
a result of significant perturbing mechanical forces, however 
the forces are not so large as to induce global unfolding. 
The mechanical forces are applied across the whole protein 
surface to ascertain a "mechanical fingerprint" for a given 
SOD1 protein variant. 

An all-atom, implicit solvent model was used, which in 
benchmark cases compared favorably to an explicit solvent 
model, but less favorably with a structure-based Go model. 
The (non-equilibrium) work values obtained by pulling a 



residue to 5A were found to strongly correlate (r = 0.96) 
with the equilibrium free energy change for the same pro- 
cess, as calculated using the weighted histogram analysis 
method (WHAM). Thus the work profiles obtained accu- 
rately represented the relative thermodynamic stability of 
various regions, or the relative thermodynamic stability be- 
tween mutant and WT for the same region. 

We found that the mechanical profile of a given SOD1 
variant could be best represented through the cumulative 
distribution of mechanical work values, where the work is 
obtained by pulling various residues in the protein out to a 
given distance. The cumulative distribution of work values 
can be thought of as a generalization of native "rigidity" , 
which accounts for the fact that the rigidity can vary place 
to place, so that the collection of work values themselves 
obey a distribution that differs between SOD1 variants. For 
purposes of distinguishing cumulative work distributions of 
SOD1 variants, we saw that mechanical data collected for 48 
residues was sufficient to represent more comprehensive sets 
of data, in that the cumulative work distribution seemed to 
have converged to within « 1 kJ/mol after a sample size of 
about 40 residues. 

Mechanical fingerprinting studies of SOD1 variants with Zn- 
binding and electrostatic loops either truncated to Gly-Ala- 
Gly linkers [38], or extended by poly-Glycine insertions re- 
vealed that, although the Zn-binding and electrostatic loops 
in apo SOD1 destabilized the /5-sandwich core of the pro- 
tein in the absence of metals, evolution has responded by 
strengthening interactions in regions flanking the loops, to 
preserve the structural integrity of the core domain and pre- 
sumably prevent misfolding. Nevertheless, mechanical fin- 
gerprinting studies of a series of C-terminally truncated 
mutants, along with an analysis of equilibrium dynamic 
fluctuations and solvent exposure while varying native con- 
straints, together revealed that the apo protein is internally 
frustrated, and that this internal strain is an allosteric con- 
sequence of evolution towards high metal-binding affinity; 
release of the stress in a truncation mutant causes loss of 
metal binding function, even though all metal-binding lig- 
ands are still present [5]. Thus, the evolutionary optimiza- 
ton of SOD1 function as a metal-binding redox substrate 
competes with apo state thermo-mechanical stability, con- 
sequently increasing the susceptibility of misfolding in the 
apo state. 

Here we examine mechanical properties of mutant and WT 
SOD1, as well as the free energetic changes accompany- 
ing post-translational modifications (PTMs) such as metal 
binding and dimerization. The affinity for metals as well as 
dimer stability are determined both by weighted histogram 
analysis methods as well as direct mechanical pulling assays 
using the Jarzynski equality, and the methods are checked 
for consistency. Of the 21 ALS-associated mutants investi- 
gated here, every one showed both reduced metal and dimer 
affinity. Some mutants had lower dimer affinity in the holo 
state than WT protein in the apo, disulfide reduced state. 
PTMs such as disulfide reduction or metal depletion also 
lower dimer stability and metal affinity, e.g. Cu-depletion 
lowers the affinity for Zn and vice-versa. Zn plays a larger 
role than Cu in determining the mechanical rigidity of the 
native state. Presence of the native C57-C146 disulfide bond 
induces stress in the apo state of the protein, which is re- 
lieved either upon metal binding, or reduction of the bond. 
Thermodynamic quantities as obtained from our simula- 
tions are compared with those obtained experimentally in 
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the Discussion section. 



Results 



2. 1 The Jarzynski/ Crooks equality overestimates the free 
energy change, due to conformational distortion 



Figure 1A plots the work performed as a function of ex- 
tension, for residue 25 of two SOD1 variants, (E,Zn) (Cu- 
depleted) SOD1 and (Cu,E) (Zn-depleted) SOD1. The pro- 
tocol is described in the Methods Section 4.1. Some variants 
have softer mechanical profiles than others; for residue 25, 
(Cu,E) is softer than (E,Zn). Figure IB shows the work to 
pull the Zn ion to a given distance from its putative position 
in WT protein, along with the free energy as obtained from 
the WHAM method, as described in Methods section 4.7.2. 
Note that the free energy is not always less than the work 
for this trajectory- it is only when ensemble-averaged that 
the thermodynamic inequality F < (W) holds. Both work 
and free energy begin to converge to their asymptotic val- 
ues at distances beyond about 5 A. Figure 1C shows the 
work to monomerize the WT dimer by pulling it apart, 
along with the free energy change as obtained from the 
WHAM method. Note that the work performed is substan- 
tially larger than the free energy change. Both work and 
free energy begin to converge to their respective asymp- 
totic values after about 10 A. The methods used for dimer 
separation are described in Sections 4.1 and 4.7.2. 

To calculate the free energy of local protein unfolding, metal 
expulsion, and dimer separation, we have calculated the 
work to either pull a particular C a to 5A, pull a metal ion 
until the force drops to zero, or pull monomers in the dimer 
apart until the force drops to zero. The various assays are 
repeated multiple times, and the corresponding free energy 
changes are obtained from the work values by applying the 
Jarzynski equality [39-41]. Finite sample-size corrections are 
accounted for [42], as described in Methods section 4.8. 

The free energy change is calculated from repeated pulling 
assays measuring the non-equilibrium work value to pull a 
residue to 5A, through the Jarzynski equality: e ~ AF / kT — 
(e~ w / kT ^ . Here AF is the free energy change, W is the 
work value, kT is Boltzmann's constant times the tempera- 
ture in Kelvin, and the brackets (• • •) denote the ensemble 
average over identical pulling assays. 

The results are shown in Figure 2: Panels (a) and (b) show 
plots of the quantity — kT In ^e — W ^ kT ") N i where the above 
ensemble average is evaluated over iV trajectories (identical 
runs), as a function of the number of trajectories N. For 
sufficiently large N, the average should converge to the 
ensemble average. The plots in Figure 2 indicate that about 
14 trajectories is sufficiently large for convergence, for the 
pulling rates and extensions we considered in our study (red 
horizontal lines in Panels (a) and (b) of Figure 2. 

Finite sample size effects reduce the estimate for the free 
energy change, however the mean dissipated work in our 
pulling assays is modest, only about a kJ/mol, indicating 




Fig. 1. (Panel A): Work performed as a function of distance 
pulled. Tethers are placed at the C a atom closest to the center of 
mass of the SOD1 monomer (H46), and the C a atom of residue 
S25 of either (E,Zn) (red) or (Cu,E) (green) SOD1, in separate 
pulling assays. Inset (a) Force-extension profile of residue 25 for 
(E,Zn) (red) and (Cu,E) (green) SOD1. Ribbon representations of 
the protein are also shown; the tethering residue is shown in licorice 
rendering (in red) and the center C a as a red sphere. The initial 
equilibrated (at A , green ribbon) and final (at 5 A , blue 
ribbon) structures are aligned to each other by minimizing RMSD; 
this indicates that overall, the structural perturbation induced by 
pulling one residue to 5 A is modest, but significant enough 
to be unlikely thermodynamically based on the magnitude of the 
work values. (Panel B): Work and WHAM-derived free energy to 
move Zn a given distance from its putative position in WT SOD1. 
Values have nearly converged after 5A. Note that the free energy 
is not always less than the work for a given trajectory (red); 
the thermodynamic inequality (AVl^) > A.F only holds true after 
ensemble averaging, as seen for the blue curve representing the 
ensemble-averaged work (over 20 trajectories) to extend the metal 
to a given distance. (Panel C): Work (red), average work (blue), 
and WHAM-derived free energy (green) to separate the WT SOD1 
dimer by a given separation distance. Values have nearly converged 
after 10A. The free energy is substantially less than the work, much 
of which goes into conformational distortion of the protein when 
harmonic constraints on C a atoms are not present. 
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Number of Trajectories Jarzynski(kJ/mol) 



Fig. 2. Comparison of free energy values calculated from Jarzyn- 
ski's equality, and from the weighted histogram analysis method 
(WHAM). (Panel a) Free energy as calculated using Jarzynski's 
equality, as a function of the number of repeated pulling assays or 
trajectories, for residue 10 to be pulled to 5A from its equilibrium 
position in a direction radially outward from the center of mass. 
(Panel b) The same plot as panel (a) for reside 17. In panels (a) 
and (b), red horizontal line indicates the value of free energy cal- 
culated from the Jarzynski Equation (2). The green line indicates 
the free energy obtained by accounting for extra corrections aris- 
ing from finite sample sizes [42] described in Method section 4.8. 
The blue line indicates the free energy value calculated by WHAM, 
wherein equilibration with umbrella sampling constraints allows con- 
formational distortion to relax, see Methods section 4.7.1. (Panel c) 
Scatter plot of the free energy values for Cu removal, as calculated 
by the Jarzynski equation (2) and WHAM with C a constraints, 
for the mutants give in Table 1. WT protein in the upper right 
of the plot is indicated, and has the largest binding free energy. 
The free energy values strongly correlate between the two meth- 
ods. (Panel d) Same as Panel (c) for the removal of Zn. (Panel e) 
Same as Panel (d), for the free energy of dimer separation. Note 
that although the free energies obtained by two methods strongly 
correlate, the slope of the best fit line is much less than unity, in- 
dicating the presence of substantial additional free energy changes 
due to conformational distortion using the Jarzynski method. (Inset 
to panel a) WHAM-derived free energy changes for metal expul- 
sion and dimer separation for WT holo SOD1, as a function of the 
stiffness of the harmonic constraints on the C a atoms. The right- 
most points on each curve correspond to the WHAM free energies 
in panels (c-e), i.e. at the values of harmonic constraint used to 
generate the WHAM free energies for those plots. The WHAM free 
energies at zero stiffness are consistent with the Jarzynski values 
of the free energy in panels (c-e). 



Table 1 

Proteins considered for mechanical scan, metal- 
affinity, and dimer-stability estimation 



Mechanical scan a /Metal expulsion/Monomerization 



SOD1 variant 


PDB ID used for 




structure generation 


WT holo 


1HL5, 2C9V 


(Cu,E) 


2R27 b ' c 


yKy U-olllll,H/ J 


2R27 b ' c 


(E Zn) 


1HL4 


dpu^on j 


1 HT ^ 1 RK7 C 


apo^ 


IrLiV < 


noio^ori J 




G127X 


1HL5 6 




1 tt'ym 


1J1Z4 V 


QUOD 


T^kl OKU 


1.T 1 V 


JJYD Y 


1 xjt 
IriL/O 


nnn A 


1 UT ce 
IrlL/O 




1 A r 7\T 


G41D 


1HL5 6 


G41S 


1HL5 6 


G85R 


2ZKW 


G93A 


3GZO 


G93C 


1HL5 6 


H43R 


1PTZ 


H46R 


2NNX C 


H46R/H48Q 


2NNX 


H80R 


3QQD 


I113T 


1UXL 


L144F 


1HL5 6 


L38V 


2WZ5 


S134N 


IOZU 


T54R 


3ECW 


W32S 


1HL5 6 







a Mechanical scans were performed for SOD1 variants above 
the horizontal line located between holo(SH) SOD1 and 
G127X. 

b Missing regions in the PDB structure were remodeled as 
described in section 4.5. 

c Structures reported have mutations from the WT sequence; 

these were "back- mutated" to construct structures for the 

WT sequence, as described in Methods sections 4.5,4.6. 

d Metal expulsion is not relevant for these variants. 

e These mutants have no PDB structure reported, so are 

constructed by mutating the appropriate residues in PDB 

1HL5. 



that pulling is near equilibrium, and thus finite sample-size 
corrections are not large: about 1-3 kJ/mol (green horizontal 
lines in Panels (a) and (b) give the refined estimate for 
the free energy change accounting for finite sample-size 
corrections). However, the free energy changes obtained from 
this method are still larger than those of the WHAM method 
(blue horizontal lines in Panels (a) and (b) of Figure 2). The 
discrepancy is about 0.6 kJ/mol when residue 10 is pulled 
to 5 A, and about 5 kJ/mol when residue 17 is pulled to 
5 A. The reason is that pulling on a particular residue at 
a finite rate distorts the rest of the protein, and there is 
a free energy change corresponding to this distortion that 
contributes to the total. The distortion may be avoided by 
applying constraints, but it is not obvious a priori how or 
where to apply such constraints when pulling residues away 
from the protein. The WHAM method minimizes the free 
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energy change due to distortion by equilibrating at each 
distance "window" in the corresponding umbrella potential. 
Convergence was tested by varying the number of windows 
between 25 and 40, and varying the equilibration time within 
each window between 10ns and 25ns, for which the results 
did not noticeably change. 

The free energy changes for metal expulsion and dimer sep- 
aration may also be calculated from both the Jarzynski 
equality and WHAM method. For these assays, distortion 
of the protein may be minimized by applying harmonic 
constraints to the Co, atoms of the protein. The harmonic 
constraints are applied in the WHAM analysis, and the re- 
sults compared with those of the (unconstrained) Jarzynski 
analysis. Side chains and other backbone atoms are allowed 
to move. 

Perhaps the main concern is whether applying such con- 
straints would impede the removal of the metal ion. The 
most direct, solvent-accessible pathway for metal expulsion 
was taken, as described in Methods section 4.2. To confirm 
that metal expulsion was not impeded by C a atom con- 
straints, the stiffness of the C a potential was varied, from 
zero (unconstrained) to 392 x 10 3 kJ/mol/nm 2 , the stiff- 
ness of a single carbon-carbon bond, as used in the WHAM 
method. If constraints impeded the process of metal expul- 
sion, a minimum would be seen as a function of stiffness of 
the constraining potential. If constraints did not impede the 
process, the free energy vs. the stiffness of the constraining 
potential would be monotonically decreasing. 

The WHAM free energy of metal expulsion and dimer 
separation for the WT protein as a function of the stiffness 
of the constraining potential is plotted in the inset to Panel 
(a) in Figure 2. The free energy is indeed monotonically 
decreasing, indicating that the constraints do not impede 
the process of metal expulsion or dimer separation. The 
free energy of deformation may thus safely be minimized by 
applying such constraints, to obtain the direct free energy 
cost of metal expulsion or dimer separation. 

Figure 2, Panels (c-e), plots the free energy changes for Cu 
expulsion, Zn expulsion, and dimer separation for the mu- 
tants of SOD1 given in Table 1. WHAM values are plotted 
on the ordinate; Jarzynski values, as given by equation (2) 
for 25 replicas (i.e. without any finite size-corrections), 
are plotted on the abscissa. The unconstrained Jarzynski 
method overestimates the free energy changes, for the rea- 
sons above. The slopes of the best fit lines are all less than 
unity, and for dimer separation the slope is particularly 
small. However the correlations between the two separate 
methods are excellent, indicating that the deformation free 
energy is not a large and random component of the to- 
tal free energy, but likely correlates with the direct metal 
expulsion or monomerization free energy. 

For both Zn and Cu binding, the ratio of the binding free 
energies, Jarzynski to WHAM, is about 1.8. Relative rank 
ordering of the free energies between the WT sequence and 
mutants may be predicted either by the WHAM or the 
Jarzynski method. We found, however, that the WHAM 
method was the most rapid and reliable way of determining 
free energy changes in the system. 

The right-most data point in Panels (c-e) refers to the WT 
protein (labelled), i.e. the WT protein has the largest metal 



and dimer affinity. This is described in further detail in Sec- 
tions 2.2, 2.3 below. Comparing the WHAM and Jarzynski 
values for say Cu expulsion in Panel (c) of Figure 2, the 
WHAM value of 55 kJ/mol is taken from the value in the 
inset to panel (a) at C a constraint of 392 x 10 3 kJ/mol/nm 2 , 
while the Jarzynski value of 100 kJ/mol is consistent with 
the unconstrained WHAM value of 99 kJ/mol in the in- 
set of Panel (a) (accounting for finite-size corrections to 
the Jarzynski estimate reduces the value by only about 0.3 
kJ/mol). Thus, the WHAM values approach the Jarzynski 
values as the C a constraints are removed, and provide a 
satisfying consistency check. 



2.2 A LS- associated SOD1 mutants show decreased affin- 
ity for both Cu and Zn 



Copper and Zinc ions were pulled out of their putative 
binding sites for the list of proteins given in Table 1. For 
mutant proteins with no PDB structure, the initial struc- 
ture was generated by mutating the appropriate residue (s) 
in the WT structure (PDB 1HL5) and equilibrating, ac- 
cording to Methods section 4.6. Some mutants of SOD1 
are devoid of Cu in the PDB structure- these are H46R, 
H46R/H48Q, S134N, D124V, H80R, and D125H, and some 
are fully metal depleted: G127X and T54R [43,44,31]. For 
these proteins, metals were incorporated into the proteins by 
superposing them onto holo SOD1. The metals were found 
to be metastable in their binding sites, i.e. no forces were 
required to constrain the metals to their putative bound 
positions for the duration of the simulation. 

The tethering residues are chosen according to Methods 
section 4.2, to select an easy pathway for metal expulsion; 
the other tethering point is the metal itself. From these 
pulling assays, 25 configurations were taken at separations 
between and 5 A, and used as initial states in umbrella 
sampling, to obtain the free energy to extract the metal 
using WHAM. 

The free energy values of metal expulsion thus obtained 
are plotted in Figure 3 (a), and also given in Table SI. 
The values in Figure 3 (a) are sorted by decreasing Cu 
affinity, and the values in Table SI are sorted by decreasing 
Cu affinity separately amongst mutants and WT PTM 
variants. The WT structure had the both the highest Cu 
and highest Zn affinity, with values of 53 kJ/mol and 40 
kJ/mol respectively. 

The free energy obtained by the WHAM method accounts 
for solvation free energy (at the accuracy of the implicit sol- 
vent model), but not fully for the chemical potential of the 
unbound metal, which is a concentration-dependent quan- 
tity. The chemical potential must be taken into account in 
determining the probability a metal is bound. Nevertheless, 
it is clear that every mutation, whether the corresponding 
amino acid was near the binding site or not, decreased the 
affinity of the protein for the metals. A non-ALS associ- 
ated mutant, W32S, also reduced the affinity of the met- 
als. The reduction in binding free energy from WT SOD1 
is up to about 26 kJ/mol for Cu (T54R) and 15 kJ/mol 
for Zn (G127X). These reductions are due to the induced 
conformational strain at the binding site that arises from 
the resulting stresses of mismatching interactions after a 
mutation. 
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. 3. (Panel a) Free energy of copper (Cu) and zinc (Zn) 
binding/expulsion for WT SOD1 and mutants, sorted by decreasing 
Cu affinity. Binding free energies are obtained using WHAM as 
described in Methods section 4.7.2. Blue bars indicate Cu binding 
free energy, red bars indicate Zn binding free energy, and green bars 
indicate monovalent Cu 1 "*" binding free energy. The heights of all 
bars start from zero. This binding energy accounts for solvation free 
energy at the level of the implicit solvent model. Inset images are 
intended to depict the direction of the pulling pathway, which was 
chosen as the direction with the most immediate solvent exposure. 
Metals are shown as an orange sphere for Cu and grey sphere for 
Zn, and the tethering residue used in the pulling simulations is 
color coded directly behind it. For Cu expulsion this was residue 
45 (blue cartoon), and for Zn expulsion this was residue 83 (red 
cartoon). (Panel b) Free energy to monomerize a homodimer of 
various SOD1 species. Both ALS-associated mutants and variants 
lacking post-translational modification (metal depleted or disulfide 
reduced) are considered. Free energies are rank ordered strongest 
to weakest: holo SOD1 has the strongest binding free energy, and 
G127X has the weakest. All mutants and variants showed weaker 
dimer stability than holo WT SOD1. The free energy to separate 
the dimer into monomers is obtained using WHAM. The inset shows 
a ribbon representation of the monomers constituting the dimer. 
The tethering residues, C a (46) on each monomer, are shown as 
red spheres, and the direction of pulling used to generate initial 
conditions in WHAM is indicated by arrows. All mutants are taken 
in the holo form, with Cu and Zn in their putative binding positions. 
As well, H46R and H46R/H48Q do not bind Cu, H80R does not 
bind Zn, and G127X does not bind either metal; in these cases we 
also calculated the binding free energy for the appropriate metal 
depleted forms. From Panel a, the Cu affinity for the T54R variant 
was the lowest of all mutants, and the Zn affinity was the 2nd 
lowest after G127X- lower than H80R which does not bind Zn. For 
these reasons, we also calculated the free energy to monomerize 
apo T54R dimer. 




Fig. 4. Schematic representations of the thermodynamic cycles 
associated with metal extraction and dimer separation of SOD1 
variants. Panel (a): Free energy changes associated with different 
steps of the metal-expulsion/metal-insertion thermodynamic cycle. 
Cu and Zn ions are shown as yellow and red spheres. AF£ ons tr 
is the free energy change for metal expulsion from the protein 
with the Cq, constraints present (see text), /s^prelax ^ g ^ e f ree 
energy change associated with the relaxation process where the 
constraints were gradually removed after metal-expulsion, /\pconstr 
is the free energy change of inserting the same metal back into the 
half-metallated protein, with C a constraints present to preserve its 
initial backbone conformation, and /^p^elax ^ g ^ e relaxation free 
energy after metal insertion. The sum of the free energy changes 
around the thermodynamic cycle is zero within the errors of the 
simulation, validating the calculation of expulsion free energies. 
Panel (b): Same as panel (a), only here the free energy changes refer 
to monomerization or dimerization compared to metal-expulsion or 
metal-insertion in panel (a). 



The free energy values of Cu and Zn expulsion are highly 
correlated, i.e. if the affinity of Cu was reduced, so was 
that of Zn. Three mutants stand out as exceptions: H46R, 
H46R/H48Q, and D125H. Residues 46 and 48 coordinate 
the Cu ion, so the Cu affinity is anomalously low in the 
corresponding mutants. The correlation coefficient between 
Cu 2+ and Zn 2 + binding free energies without these 2 mu- 
tants was 0.929, with significance P r = 4.4e-10. The correla- 
tion coefficient between monovalent Cu 1+ and Zn 2 + binding 
free energies without these 2 mutants was r = 0.930, with 
significance P r = 3.7e-10. The mean occupation numbers of 
Cu and Zn for several SOD1 mutants have been obtained 
experimentally by Hayward et al [45]. Re-analyzing this data 
to inspect the correlation between Cu and Zn content gives 
a correlation of r = 0.87 (p = 5e-6), r = 0.84 (p = le-4) if 
mutants H46R, H48Q, and D125H are removed. The strong 
effect of D125H on Cu affinity is more subtle and involves 
propagation of strain through the native protein; further 
description of the effect is given in the Discussion section. 

SOD1 catalyzes the disproportionation of superoxide anion 
(O^ - ) in a two-step reaction, wherein the Cu cation cycles 
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between divalent and monovalent states. In the monovalent 
state, the binding affinity for Cu + is reduced by slightly 
more than a factor of 2, resulting in more likely release 
of the cation, particularly for mutants. In fact, affinities 
for monovalent Cu+ were always less than those of Zn 2 +, 
indicating ready release of monovalent Cu upon structural 
perturbation or rare thermodynamic fluctuation. 

We have further checked that metal depletion and monomer- 
ization each satisfy a thermodynamic cycle, such that 
upon metal re-insertion or re-dimerization, along with re- 
equilibration after the respective process, the net free energy 
change over the cycle is zero within the errors of the simu- 
lations. The thermodynamic cycles are shown schematically 
in Figure 4, and the free energy changes for the respective 
processes are tabulated in Table S2 of the Supplemen- 
tary Content. On average ^2 cycle AF t « 0.15 kJ/mol for 

metal expulsion/insertion, and f AFi « 0.2 kJ/mol 

for monomerization/dimerization, for the cases that we 
benchmarked. 



2.3 ALS-associated SOD1 mutants show increased ten- 
dency to monomerize 



Dissociation of the native SOD1 dimer is a prerequisite for 
its aggregation [10,46,47]. The pathway to aggregation gen- 
erally proceeds by non-native dimer formation and subse- 
quent oligomer formation. We investigated how both mu- 
tational and post-translational modification modified dimer 
stability, for the ALS-associated mutants and PTM variants 
of SOD1 given in Table 1. 

Using the WHAM methodology described in Methods sec- 
tion 4.7.2, we found the free energy to separate the ho- 
modimer into a pair of monomers. A plot of the dimer 
binding free energies, rank ordered from strongest to weak- 
est, is shown in Figure 3 (b). The absence of any post- 
translational modification (e.g. absence of a metal, reduction 
of the disulfide bond) lowered the free energetic stability 
of the dimer. All ALS-associated holo mutants had reduced 
dimer stability, and several that we investigated had less 
dimer stability than that of apo(SH) WT, which is known 
to monomerize [11]. On the other hand, apo(SH) WT SOD1 
had lower dimer stability than 15/22 holo ALS-associated 
mutants. The non- ALS-associated mutant W32S showed re- 
duced dimer stability. G127X, an obligate monomer in in 
vitro studies [31], had the lowest dimer stability: a signif- 
icant fraction of the residues participating in the dimer 
interface is removed by the terminal sequence mutation. 

The trend in dimer stability did not correlate with the prox- 
imity of a mutated residue to the dimer binding interface, 
either by measuring the distance to the closest residue in 
the binding interface (r = 0.15,P = 0.51), or the mean dis- 
tance to residues putatively involved in the dimer interface: 
residues 4, 50-53, 114, 148, 150-153 (r = 0.096,P = 0.67). 
These results support a dimer destabilizing mechanism by 
mutations that involves the long-range propagation of stress 
through interaction networks in the protein, an interpreta- 
tion that is consistent with previous experimental [48] and 
simulation studies [49]. 




Work(kJ/mol) 



Fig. 5. The effect of metal depletion on the mechanical stability 
of SOD1. (Inset panel a) Mechanical profiles of holo SOD1 (black), 
(E,Zn) SOD1 (blue), (Cu,E) SOD1 (green), and (Cu,E) SOD1 with 
the Cu ion shifted from its putative position to that of the Zn 
(red), which is the location of the metal observed in the crystal 
structure of PDB 2R27 [7]. Profiles are color-coded according to the 
legend in the upper left. The effect of metal depletion modulates 
the entire mechanical stability profile. Work values for individual 
residues are given in Table S5 of the Supplementary Content. (In- 
set b) Probability distributions of mechanical stability for holo and 
(Cu,E) variants. This representation shows that the distributions are 
indeed different, however there is significant overlap, which along 
with the necessity of binning the data to construct the histogram, 
makes the differences difficult to quantify. (Main Panel) Cumulative 
distributions of the above SOD1 variants. The cumulative distribu- 
tion requires no binning, and simply rank orders and accumulates 
the work values. Using this representation, discrepancies between 
work profiles most clearly emerge. The mechanical stability of the 
SOD1 variants can be rank ordered, strongest to weakest, as holo, 
(E,Zn), (Cu-shift,E) , and (Cu,E). Also shown in the main panel is 
a correlation matrix resulting from a pairwise comparison of the 
work values for all pairings of SOD1 variants. The profiles do not 
correlate, again supporting the notion that metal depletion mod- 
ulates the entire mechanical stability profile. A similar plot also 
containing the cumulative distribution for apo SOD1 is given in 
the Supplementary Content (Figure SI). 

2.^- Metal depletion weakens the mechanical stability of 
WT SOD1- Zinc moreso than Copper 

To investigate the effect of the presence or absence of Cu 
and Zn ions on the mechanical stability of WT SOD1, 
we examined the mechanical profiles of the metal present 
(holo), Cu-depleted (E,Zn), Zn-depleted (Cu,E), Zn-depleted 
with Cu shifted to the Zn position (Cu-shift,E), and metal- 
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depleted (apo) forms of SOD1. Mechanical profiles were 
obtained from pulling simulations as described in Methods 
section 4.1. Inset (a) of Figure 5 shows the mechanical 
profiles of all the above variants except for apo SOD1 (an 
analogous figure including the data for apo SOD1 is given 
in Figure SI in the Supplementary Content). 

The work values in the mechanical profile are significant 
compared to thermal energies, ranging from about 16/c#T 
to Though the structural perturbations are fairly 

modest in comparison to global unfolding, the energetics of 
the pulling process is sufficiently non-perturbative that the 
work values across the set of 48 residues used here do not 
correlate with the RMSD values for the respective residues 
(see Table S3 in the Supplementary Content). 

The mechanical profiles are significantly different for all the 
variants- a table providing the cross-correlation coefficients 
is given to the left of Figure 5, and also in Table S4 in the 
Supplementary Content, which shows that no mechanical 
profiles is well-correlated with any other. It is worth not- 
ing however that the correlation coefficients, though weak, 
are statistically significant between some variants: The p- 
value between (Cu,E) and (Cu-shift,E) is 0.015, and the 
p-value between holo and (Cu-shift,E) is 0.028. The effect 
of demetallation is difficult to discern from the mechanical 
profiles in inset (a). Here the cumulative distribution proves 
to be valuable in elucidating the mechanical discrepancies 
between SOD1 variants, and is shown in the main panel of 
Figure 5. The probability distribution of work values can 
differentiate work profiles for different SOD1 species (in- 
set (b) of Figure 5), but unlike the cumulative distribution 
it requires binning, and does not discriminate variants as 
clearly as the cumulative distribution. 

According to the cumulative distributions in Figure 5, metal 
depletion mechanically destabilizes the protein. That is, af- 
ter depletion of either Cu or Zn, the protein is more me- 
chanically susceptible to perturbing forces. The distribution 
of apo SOD1 is broadened compared to (Cu,E) SOD1, but 
not overall weakened (Figure SI, Supplementary Content). 
The relative importance of Cu and Zn in determining the 
mechanical rigidity of the protein can also be ascertained. 
Zn depletion results in more destabilization than does Cu 
depletion, over the whole range of work values observed 
(green and blue curves in Figure 5). That is, at any given 
perturbing work value, more residues in the (Cu,E) pro- 
tein would be substantially disordered than in the (E,Zn) 
protein. Equivalently, for say the weakest 1/3 (or any frac- 
tion) of the residues, a smaller value of work is required to 
substantially disorder those residues in the (Cu,E) variant 
than in the (E,Zn) variant. Note that the residues in the 
weakest 1/3 may be different in the two variants. 

The loss of Cu and/or Zn destabilizes both the Zn-binding 
loop and the electrostatic loop, Zn moreso than Cu. Zn 
has close proximity to some of the residues in the elec- 
trostatic loop, and is also involved in helix dipole capping 
interactions with helix 133-138 [50]. From the work values 
in the mechanical profiles (See Table S5 in the Supplemen- 
tary Content), the mean work for residues in the mechan- 
ical scan that were in the Zn-binding loop (residues 50, 
54, 55, 60, 65, 70, 73, 75, 80) was 83.9 kJ/mol for holo, 
78.2 kJ/mol for (E,Zn), 76.4 kJ/mol for (Cu,E), and 71.5 
kJ/mol for apo. The mean work values for the correspond- 
ing residues in the electrostatic loop (residues 121, 125, 130, 
135, 136, 140, 141) was 93.5 kJ/mol for holo, 85.6 kJ/mol 



for (E,Zn), 77.7 kJ/mol for (Cu,E) and 72.9 kJ/mol for apo. 
The changes in work values with respect to holo SOD1 are 
largest for the zinc binding and electrostatic loops (Figure 
S7), indicating that the mechanical stability of these regions 
is preferentially dependent upon metal content. This find- 
ing is supported by dynamical fluctuation data, which show 
preferential increase in RMSF values for the zinc binding 
and electrostatic loops (Figure S7 and Table 2). 

The changes in the mechanical stability profiles as a function 
of sequence index, when metals are removed, are subtle to 
predict and would likely involve a detailed quantification 
of the network of interactions throughout the protein. The 
modulus of the changes in work values do not correlate with 
simple parameters such as distance to the either metal (all 
correlation coefficients had r < 0.26). 

In the crystal structure of (Cu,E) SOD1, Cu resides near the 
putative Zn position. We tested whether this shift in position 
lowered the free energy by pulling the Cu from its putative 
position in the crystal structure to the Zn binding position, 
and applying umbrella sampling with the WHAM method as 
described in Methods section 4.7.3, to obtain the free energy 
change for such a shift. Indeed these simulations gave a 
reduction in free energy of —5.7 kJ/mol, indicating that the 
Zn binding position is more favorable for the Cu ion, when 
only Cu is bound. However the free energy barrier between 
the Cu and Zn binding positions is sufficiently large, about 8 
kJ/mol, that we did not see the Cu change binding positions 
during the course of our simulations (covering 0.2 /is) when 
Zn was absent. We tested the hypothesis that the shift of 
Cu also mechanically stabilizes the protein, by calculating 
the mechanical profile and cumulative distribution of (Cu- 
shift,E) SOD1. These results are shown in Figure 5, and the 
values are given in Table S5 of the Supplementary Content. 
Shifting of the Cu ion to the Zn position significantly 
increases in the mechanical stability of the protein: The 
probability to obtain, by chance, a distribution as stable 
or more stable than the (Cu-shift,E) distribution from the 
(Cu,E) distribution is w 1.4e-7 (see Methods Section 4.9 and 
Figure S2 in the Supplementary Content). As mentioned 
above, the work profiles of (Cu,E) and (Cu-shift,E) SOD1 
are only weakly correlated: r = 0.35, p = 0.015. (correlation 
table in Figure 5). Shift of the Cu to the Zn position globally 
changes the mechanical malleability of the protein. 

The work profiles of both (Cu,E) SOD1 and (Cu-shift,E) 
SOD1 show lower stability in the metal binding residues 
H80 and H120, compared to holo SOD1. The work values 
for residue 120, which coordinates the Cu 2+ ion, are 81.6 
kJ/mol in the holo state, 56.6 kJ/mol in the (Cu,E) state, 
and 44.7 kJ/mol in the (Cu-shift,E) state. That is, Zn- 
depletion destabilizes Cu-binding residue 120, and shift of 
the Cu ion to the Zn-binding position further destabilizes 
that residue. Similarly, the work values for residue 80, which 
coordinates Zn 2 +, are 73.1 kJ/mol in the holo state, 57.6 
kJ/mol in the (Cu,E) state, and 68.6 kJ/mol in the (Cu- 
shift,E) state. That is, Zn-depletion destabilizes coordinating 
residue H80, and shift of the Cu to the Zn-position, a 
free energetically favorable process, partially recovers the 
mechanical stability of that residue in the WT protein. 

The distribution of apo SOD1 is broadened compared to 
(Cu,E) SOD1, but not weakened overall, (the mean work 
for both variants was within 0.2 kJ/mol). That is, the 
most weakly mechanically stable residues in apo SOD1 are 
less stable than the weakest mechanically stable residues 
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in (Cu,E) SOD1 (though these residues are not the same), 
but the most stable residues of apo SOD1 are more stable 
than the most stable in (Cu,E) SOD1 (and these residues 
are also not the same in the two variants). 



2.5 Metal binding and disulfide bonding are cooperative: 
removing one destabilizes the other 



Figure 6 plots the cumulative distributions of the mechanical 
work profiles for holo, apo, holo/SS-reduced (holo(SH)), and 
apo/SS-reduced (apo(SH)) SOD1. All modifications desta- 
bilize holo SOD1. By comparing apo with holo(SH) SOD1 
however, metal depletion has the largest destabilizing ef- 
fect with only a few exceptions, over the bottom 77% of 
the work values. The cumulative distributions then cross at 
around 81 kJ/mol. 

The weakest regions, with work values less than 60 kJ/mol, 
are on average about 5 kJ/mol weaker in apo SOD1 than 
in holo(SH) SOD1, while the residues with moderate me- 
chanical stability in the range of 60 — 80 kJ/mol are only 
about 1.4 kJ/mol weaker on average for apo SOD1. In this 
regime, disulfide reduction has nearly as much destabilizing 
effect as metal depletion. The most mechanically stable re- 
gions, with work values > 86 kJ/mol, are on average about 
5.5 kJ/mol weaker for holo(SH) SOD1 than for apo SOD1, 
i.e. disulfide reduction is more destabilizing for the most 
mechanically stable residues requiring the highest work val- 
ues. It is worth emphasizing again that these regions in the 
cumulative distribution need not involve the same residues 
for the 2 variants. 

Perhaps surprisingly, metal depletion appears to have a 
moderately stabilizing effect on holo(SH) SOD1, at least 
over the range of work values less than 81 kJ/mol. Apo(SH) 
SOD1 contains the 3 weakest residues between the holo(SH) 
and apo(SH) variants, however for 19 out of the total of 31 
residues of the mechanical scan with work values less than 81 
kJ/mol, apo(SH) SOD1 is more stable than holo(SH) SOD1. 
This is seen directly from the cumulative distributions. This 
observation prompted a study of the mechanical rigidity 
of SOD1 variants in the vicinity of the residues involved 
in the disulfide bond, for various metallated states. For a 
given metallated state of the protein, we recorded the mean 
mechanical work to pull residues 57 and 146 (involved in 
the disulfide bond) to 5 A, for both the disulfide-present 
and the disulfide-reduced state. If the disulfide bond is 
present, both residues will tend to move together when one 
is pulled. Inset (a) of Figure 6 plots the difference, SS- 
present minus SS-reduced, of the mean mechanical work to 
pull residues 57 and 146, for holo, (E,Zn), (Cu,E), and apo 
SOD1. For holo SOD1, the presence of disulfide bond results 
in mechanical stabilization of these residues, or stronger 
coupling to the rest of the protein. The effect goes away 
once Cu is removed from the protein. If Zn is removed, the 
reverse effect is observed: formation of the disulfide bond 
weakens the mechanical coupling between residues 57/146 
and the rest of the protein. The mechanical weakening effect 
of the disulfide bond is most significant in the apo state of 
the protein. 

Just as metal depletion mechanically stabilizes the disulfide- 
reduced state of the protein, so also does disulfide-reduction 
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Fig. 6. Interplay of disulfide linkage and metal depletion on 
protein stability. (Main panel) Cumulative distributions of work 
values for holo, apo, holo(SH), and apo(SH) SOD1. Both metal 
depletion and disulfide reduction reduce the stability of the WT 
protein, but metal depletion generally results in larger destabilization 
for work values less than about 80 kJ/mol. Comparison of the holo, 
holo(SH), and apo(SH) cumulative distributions indicates that the 
presence of the disulfide bond stabilizes the holo form of SOD1, 
but destabilizes the apo form of SOD1. (Inset a) Bar-plot of work 
differences for several SOD1 variants using the following procedure. 
We first found the average work, W , to pull residues 57 and 146 
to 5 A, for both disulfide-containing and disulfide-reduced forms 
of holo, (E,Zn), (Cu,E) and apo (metal-depleted) SOD1. Then we 
took the difference in average work between the disulfide present 
and disulfide-reduced forms for all those variants W$S ~ ^Red- 
This was taken as a measure of effect of disulfide bond on the 
stability of SOD1. The disulfide bond stabilizes the holo form, but 
destabilizes the (Cu,E) and apo forms. (Inset b) Potential energy 
as a function of time after initial reduction of the disulfide bond at 
t = 0. Red curve is for holo SOD1, green curve is for apo SOD1. 
(Inset c) Schematic free energy profiles of apo and apo(SH) SOD1 
(see Discussion). The key features of the profiles are a stiffer native 
basin, a moderately lower free energy in the native state, and a 
substantially lower free energy in the unfolded state, for apo(SH) 
SOD1 relative to apo SOD1. The mechanical profiles probe the 
region outlined in the dashed box. 

stabilize the metal-depleted protein. Comparing the cumu- 
lative distributions of apo(SS) and apo(SH) SOD1, 32 of 38 
work values less than 86 kJ/mol are more stable for apo(SH) 
than for apo(SS). Thus disulfide-reduction mechanically sta- 
bilizes at least the weaker regions of the apo protein. 

What about the effect on the overall potential energy in 
the protein? We considered the following kinetics study: in 
both holo and apo SOD1, we instantaneously reduced the 
disulfide bond, and investigated the potential energy in the 
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protein as a function of time. Inset (b) of Figure 6 plots 
the potential energy as a function of time immediately after 
disulfide reduction, for the holo protein (red), and the apo 
protein (green). The holo protein indeed shows an increase 
in potential energy indicating a destabilizing effect due 
to disulfide reduction. However, consistent with the above 
observations from the cumulative distribution and the local 
work on disulflde-participating residues, the potential energy 
in the apo protein decreases with time after reduction. This 
indicates that the presence of the disulfide bond strains the 
apo protein, and the potential energy in the protein is thus 
lowered by its reduction. 

The converse study of removing metals from the protein 
with the disulfide bond either present or absent is con- 
tained in Figure 3 (a) and Table SI. The holo(SH) protein 
has moderately reduced affinity for Cu and Zn, by about 
5.9 kJ/mol and 3.9 kJ/mol respectively. Thus, reducing the 
disulfide bond lowers the affinity for the metals, and remov- 
ing the metals changes the effect of disulfide bond formation 
from protein-stabilizing to protein-destabilizing. 



3 Discussion 



Operating on the premise that ALS-associated mutant se- 
quences of SOD1 protein had different mechanical properties 
than the WT protein, we have undertaken a computational 
study to probe these mechanical differences. By employing 
a combination of pulling simulations and umbrella sampling 
with the weighted histogram analysis method (WHAM), 
we can computationally investigate disease-relevant misfold- 
ing processes such as loss of native structure, metal loss, 
and monomerization of native dimers. We found the likeli- 
hoods of these processes to be significantly greater in ALS- 
associated mutant proteins, supporting the idea that al- 
though the consequences for folding thermodynamic stabil- 
ity of many ALS-associated mutants may be modest [51,45], 
the effects on dimer stability and metal affinity can be 
large. The commonly-used monomeric F50E/G51E mutant 
provides a clear precedent for such effects. Moreover, dimer 
and metal affinities can be influenced in subtle ways: mu- 
tations distant from the Zn or Cu can have large effects on 
metal affinity, and mutations far from the dimer interface 
can have large effects on dimer stability. 

Mechanical probes were employed by simulating tethers on 
the center residue of SOD1, and on various residues on 
the protein surface. The experimental analogue to such an 
in silico approach would require multiple AFM or optical 
trap assays involving numerous residue pairs about the 
protein surface as tethering points. This is difficult and 
time-consuming to achieve in practice, which thus provides 
an opportunity for the present simulation approaches. Such 
a study of the native surface malleability for mutants of 
SOD1 may be relevant to understanding the process of 
intermolecular transmission of misfolding in the cell [31]. 

The work to pull a given residue to a distance sufficient to 
constitute an anomalously large fluctuation, e.g. 5A, can be 
calculated as a function of sequence index for a given SOD1 
variant. This results in a characteristic mechanical profile for 
that protein. We found that such a profile was significantly 
different between WT SOD1 and other post-translationally 



modified (PTM) variants such as metal depleted or disul- 
fide reduced SOD1. A systematic comparison of mechanical 
profiles between WT SOD1 and ALS-associated mutants is 
an interesting and important future study. 

Some of our in silico observations are consistent wth 
previous experimental benchmarks, however we have also 
made several new observations here that are experimentally 
testable. Table 2 gives a summary of the experimentally- 
validated results in this paper, as well as the main 
experimentally-testable predictions contained in the paper. 
We elaborate on the entries in this table in the discussion 
below. 

A given PTM globally modulates the mechanical profile, in- 
ducing both local and non-local changes in stability. These 
changes can be destabilizing in some regions and stabilizing 
in others. Understanding such stability changes well enough 
to predict them is a difficult challenge, and would involve 
quantifying the interaction networks throughout SOD1, and 
how these networks transmitted stresses when one region 
of the protein was strained. A detailed study of the con- 
sequences of such mutations, for example by studying the 
long-range communication through interaction networks in 
the mutant vs WT protein (see for example Khare and 
Dokholyan [49]) is an interesting topic of future research. 

The free energy to remove metals or monomerize SOD1 vari- 
ants may be calculated using WHAM. For these processes, 
harmonic distance constraints between C a atoms were ap- 
plied to prevent the inevitable structural deformations that 
occur in such a non-equilibrium pulling assay, which would 
add to the measured free energy. The free energy as ob- 
tained from the Jarzynski equality (without distance con- 
straints) was thus always greater than the free energy ob- 
tained from the WHAM protocol along with harmonic dis- 
tance constraints. Removing the distance constraints from 
the WHAM protocol reproduced the Jarzynski-derived free 
energy. 

We found that many discrepancies of SOD1 variants from 
WT were difficult to disentangle from the work profile, but 
often emerged most naturally from the cumulative distri- 
bution of work values. Mechanical scans were thus used to 
construct the cumulative distributions, which then allowed 
us to distinguish the stabilizing energetics in various forms 
of PTM SOD1. In principle, different mechanical profiles 
could give rise to similar cumulative distributions, but in 
practice this was never an issue. If cumulative distributions 
for two variants happened to be similar, their mechanical 
profiles could always be compared. Many of the proper- 
ties we may be interested in for a given SOD1 variant in- 
volve coarse-grained properties of the protein that can be 
directly obtained from the cumulative distribution, such as 
how many residues are weakly stable given a threshold of 
work that perturbs the system. 

By comparing the cumulative distribution of the work values 
for holo and metal depleted SOD1 variants, we found that 
loss of either metal makes the protein more susceptible to 
structural deformation. Loss of Zn had a larger effect than 
loss of Cu. The loss of either metal modulated the entire 
mechanical profile, such that by pairwise comparison, the 
work values no longer showed strong correlation with those 
in the holo protein. The strongest correlation with holo 
protein was seen for (Cu-shift,E) SOD1, which lacks Zn and 
has Cu shifted to the Zn position (r = 0.32, p = 0.03). 
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Table 2 

Experimental validation and predictions of the model 



Property 


Simulation 


Experimental Reference 


AG of apo WT dimer — y apo monomer 


62 kJ/mol (Fig. 3b) 


61 kJ/mol [52] a 


ACr of apo G85R dimer — y apo monomer 


4/.y KJ/moi \r lg. ODJ 


A Q o ~\r T /vy-.^1 \AQ~\ 

4o.y KJ/moi [4oj 


Non-local destabilization of 
dimer stability of SOD1 mutants 


Fig. 3b (and text) 


[48], [49] b 


Zn-binding AG for WT SOD1 


40 kJ/mol (Fig. 3a) 


/in \r T /rvi/~J T/l^l /I £ 1^ T /w-i^l T^QlC 

4U KJ/moi [4DJ , 40 KJ/moi [OoJ , 
56 kJ/mol [12], 76 kJ/mol [54] 


I-. i I-. a TA^ r r oath 
LyU-Dinamg ZaLt tor w l oiJJJl 


oo KJ/mol \r lg. oaj 


yo KJ /mol [04J 


Zn-binding A AG for mutant SOD1 
Rvalue (^mutant ) ) r — u.yy,p — u.ui 


3.2(A4V), 4.6(G93A), 

( . ( {Lioov ) , 1Z.0(^o1041N ) KJ/mol 


-2.9(A4V), -1.7(G93A), 

Z.y^L/ooVJ, 0./^olo41NJ KJ/mol [IzJ 


Cu-binding AAG for mutant SOD1 
(value(mutant)) r=0.78 d 


8.1(A4V), 14.9(I113T), 
12.0(138V) kJ/mol 


1.2(A4V), 2.1(I113T), 
2.4(L38V) kJ/mol [54] 


uu,rj bODl is more malleable than E,Zn bODl 
(Zn-binding facilitates native 
structure moreso than Cu binding) 


Fig. 5 


[36,55-58] 


Zn-binding and electrostatic loops 

are selectively destabilized by metal release 


Fig. S7 


[7,13,14,36,58] 


apo(SS) is less stable than holo(SH) SOD1 


Fig. 6, Fig. S6, Table S5 


[59] 


"\ T 1*1 j • £ T 1 * 1 • i • i_ 1 

Validation 01 Jarzynski equality with 
WHAM for native mechanical stability 


Fig. 2 


[42] 


Predicted Property 


Simulation 


Suggested Experiment 


Dimer stability tor bUUl variants and ALb 
mutants e (all had reduced AG w.r.t. WT; 
D124V, H80R, and G85R were among lowest) 


Fig. 3b 


Equilibrium denaturation 


Zn-binding/Cu-binding AG for SOD1 variants 
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a Dimer binding free energies for pseudo-WT C6A/C111S are somewhat less, about 51-52kJ/mol/dimer [52,48,61]. 
b Simulation result. 

c This number is reported as an approximate lower bound to the affinity. 
d Dataset is too small (N=3) for the correlation to have statistical significance. 

e Mutants with newly predicted dimer binding free energies: holo(SH), (Cu-shift,E)(SS), Cu,E(SS), E,Zn(SS), E,Zn(SH), the 
following apo mutants: T54R, H80R, and G127X, and the following holo mutants: A4V, G37R, L38V, G41D, G41S, H43R, 
H46R, H46R/H48Q, T54R, D76Y, H80R, G93C, I113T, D124V, D125H, G127X, S134N, L144F. 

Mutants and variants with predicted metal binding free energies for Zn: holo(SH), E,Zn(SS), E,Zn(SH), W32S, G37R, 
G41D, G41S, H43R, H46R, H46R/H48Q, T54R, D76Y, H80R, G85R, G93C, D124V, D125H, G127X, L144F; for Cu: 
holo(SH), (Cu-shift,E)(SS), Cu,E(SS), W32S, G37R, G41D, G41S, H43R, H46R, H46R/H48Q, T54R, D76Y, H80R, G85R, 
G93A, G93C, D124V, D125H, G127X, S134N, L144F. 

^Proteins used: WT, A4V, L38V, G41S, G72S, D76Y, D90A, G93A, E133A, H46R, H48Q, G85R, D124V, D125H, S134N 
(numbers in parentheses do not include H46R, H48Q and D125H). 

^The fact that NMR measurements report no significant change in tertiary structure upon disulfide reduction of monomeric 
apo(SS) SOD1 [8] is consistent with this prediction. 
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Interestingly, the work values of apo SOD1 showed mod- 
est anticorrelation with those in holo SOD1 (r = —0.37, 
p = 0.01), consistent with the idea that mechanically rigid 
residues in the holo protein are stabilized by the presence 
of metals, and that mechanically rigid residues in the apo 
protein tend to be destabilized by the metals. This compe- 
tition between apo-state stability and metal affinity-driven 
holo-state stability was also seen in mutant studies, which 
revealed significant strain in the apo state [5]. 

Our simulations showed that (Cu,E) SOD1 is substantially 
more mechanically malleable than (E,Zn) SOD1, which is 
consistent with previous experimental data: Early work by 
Rotilio et al found that the extent of recovery of the na- 
tive copper site in terms of spectroscopic features and cat- 
alytic activity was related to the amount of residual zinc 
present [55], implying that Zn aided the incorporation of 
Cu. Work from the Valentine group later showed, by an- 
alyzing NMR shifts of coordinating histidine residues in 
the presence of a chemical modifier contingent on solvent 
exposure [56], that Zn aided the incorporation of Cu by 
structuring the protein, and that Zn alone was sufficient 
to structure the protein. Matthews and co-workers found 
substantial stabilizing free energies upon Zn binding, for 
all phases of structural organization of SOD1, dimeric to 
unfolded [57]. 

Crystal structures of mutants that either inadequately or 
anomalously bind Zn show amyloid filament-like non-native 
interactions, further supporting the critical structural role 
of Zn. Elam et al [36] find that the crystal structures of apo 
H46R, and S134N, for which Chain B has only 1 Zn ion 
bound in the Cu position, both show similar and significant 
structural disorder in which non-native contacts are made 
between the electrostatic loop from one domain and an 
exposed cleft made by {3 strands 5 and 6 of a neighboring 
domain. Oliveberg and co-workers find a similar effect for 
the H63S/H71S/H80S/D83S variant [58], which cannot bind 
Zn in the putative binding site, but instead binds Zn in the 
Cu binding site for about 3/4 of the molecules. Histidine 
48 fails to ligate the Zn ion and instead swings outwards, 
leading to structural disorder in loop IV, which in turn 
destabilizes Loop VII. This loop then interacts with edge 
f$ strands in neighboring molecules, analogously to the apo 
H46R and S134N crystal structures. 

Consistent with the observations that the role of Zn is struc- 
tural while the role of Cu is enzymatic, we found that the 
affinity for Cu is reduced by ~ 11.4 kJ/mol in the absence of 
Zn, while the affinity for Zn is reduced to a lesser degree by 
« 5.9 kJ/mol in the absence of Cu. In the simulation pro- 
tocol we employed, all structures were equilibrated before 
applying constraints between Ca atoms, so the structural 
deformations associated with loss of the partner metal prior 
to the affinity measurement were accounted for. As well, 
we also accounted for the effects due to re-equilibration of 
the protein structure in the final state of the protein af- 
ter extraction of the metal, in calculating the net affinity. 
This relaxation effect is more significant for Zn than for 
Cu due to the stronger coupling of native structure to Zn 
binding, and most significant for dimer separation, because 
of the removal of steric constraints and large change in in- 
teraction energies (see Figures S3, S4, and S5). However, 
metal binding free energies were obtained from a metal ex- 
pulsion process along "the most likely" pathway, but they 
did not consider all possible pathways of escape- the con- 
tribution of all other pathways is an entropy contribution 



that would likely result in modest reduction in the estimate 
for the binding free energy. Nevertheless, our measurements 
should give a good approximation to the relative differences 
in affinity between species. The decreased metal affinity of 
ALS-associated mutants points towards the general role of 
metal loss in the misfolding and propagation process. 

Matthews and co-workers obtain larger values of Zn binding 
free energies of about 56 kj/mol/monomer [12]. Although 
these values are obtained from several approximations, in- 
cluding the use of C6A/C111S and C6A/C111S/F50E/G51E 
as reference states for dimeric and monomeric WT SOD1, 
it is instructive to investigate reasons on the computational 
side as to why their value might be larger than our value of 
40 kJ/mol. Our calculations obtain the binding free energy 
for monomeric SOD1, and so neglect cooperativity effects 
between monomers in Zn binding to the dimer state. As 
well, although the OPLS-AA/L force field used in our sim- 
ulations attempts to account for the electronic polarization 
of Histidines which are critical in modulating metal binding 
free energies, polarization may still be underestimated by 
the force field. The discrepancy between our value of Cu- 
binding free energy to WT SOD1 and that of Beckman and 
colleagues [54] indicates that Histidine polarization may be 
particularly important in Cu-binding. 

Further information on the loss of structure due to metal 
depletion is obtained from studies of the solution structure 
of the obligate monomer mutant E133Q/F50E/G51E, for 
which Banci et al. [13] found that the electrostatic and zinc- 
binding loops were severely disordered and extensively mo- 
bile in the absence of both metals, but not in the absence of 
Cu alone, i.e. for the (E,Zn) variant. Crystal structural stud- 
ies of the apo WT protein also confirm that for monomers 
devoid of metals, the electrostatic and Zn-binding loops are 
disordered and not visible in the electron density maps [14]. 
Some monomers in this study with 20% occupancy at the 
Zn site possessed well-ordered Zn-binding and electrostatic 
loop regions, however this was thought to be due to crystal 
packing forces. Beckman and colleagues have crystallized a 
constitutively Zn-deficient but Cu-containing SOD mutant 
H80S/D83S/C6A/C111S by mutating two zinc-binding lig- 
ands to hydrogen-bonding serines [7]. This variant was ob- 
served to have a disrupted dimer interface and partly dis- 
ordered zinc binding and electrostatic loops. Both our me- 
chanical pulling and simulation results are consistent with 
these observations (Table 2). In Figures S7b and c, the in- 
crease in mechanical malleability and magnitude of fluctu- 
ations upon metal release are measured for several regions 
SOD1, including residues 1-48, residues 48-83 constituting 
the Zn-binding loop, residues 84-120, residues 121-142 con- 
stituting the electrostatic loop, and residues 143-153. These 
plots show that loops IV and VII undergo systematic soft- 
ening and increased fluctuations compared to other regions, 
as metals and in particular Zn are lost from the protein. 

The mechanical response to the loss of metals that we 
observed was not localized to loops IV and VII, and in this 
sense apparently differs from the above observations and 
constitutes a prediction of the model (Table 2). For example, 
while the mean mechanical stability of loops IV and VII 
decreased substantially after the loss of ZN or both metals, 
some regions of the protein such as strands 3 and 6 are 
more rigid for the apo and (Cu,E) proteins than for the holo 
protein. Metal depletion of the WT enzyme modulated the 
entire mechanical stability profile, locally and nonlocally. 
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Our pulling simulations, which incur a significant perturba- 
tion on the structure, measure a quantity that is distinct 
from thermal fluctuations in the native ensemble. As such, 
they may show discrepancies with experimental probes mea- 
suring equilibrium fluctuations. Supporting this, the work 
values to extend to 5A do not correlate with the RMSD 
values for the corresponding residues (Table S3, Supple- 
mentary Content). As well, normal modes of the fluctu- 
ation spectrum for beta sheets may involve, for example, 
sliding motions that are transverse to the radially-directed 
pulling simulations that we employed. On the other hand, 
misfolding-related processes likely involve motions that are 
larger deformations beyond native basin fluctuations, so the 
information gained from relatively large mechanical pertur- 
bations can be useful in understanding misfolding. 

We saw that post-translational modifications, as well as 
ALS-associated mutations, reduced both metal affinity and 
dimer stability. This is consistent with experimental ob- 
servations that mutations seemed to favor increased for- 
mation of Zn-free monomeric intermediates [62]. Several 
ALS-associated mutations that we investigated reduced holo 
dimer stability moreso than did disulfide reduction or loss 
of metals from the holo WT protein. On the other hand, 
apo(SH) WT SOD1 had lower dimer stability than 15/22 
holo ALS-associated mutants; metal depletion along with 
disulfide reduction has larger consequence for monomeriza- 
tion than many mutations do alone. The apo mutants may 
have lower dimer stability still, however. The decrease in 
dimer stability of holo WT protein due to disulfide reduc- 
tion, about 37 kJ/mol, is more significant than depletion 
of either metal, but not both (about 43 kJ/mol). The de- 
crease in dimer stability due to Cu depletion is the smallest 
change of all modifications we investigated, about 12 kJ/mol, 
while the decrease in dimer stability due to Zn depletion 
is significant: about 28 kJ/mol. Increased concentrations of 
oxidatively-modifled monomeric species have been observed 
on pathway to aggregate formation by dynamic light scatter- 
ing measurements [10], supporting the notion of monomers 
as an amyloid precursor. As well, aggregation rates have 
been observed to increase with decreased concentration of 
SOD1, which strongly indicates that monomeric species act 
as precursors to aggregation [46]. Small molecules targeted 
to stabilize the dimer interface have been observed to signif- 
icantly inhibit in vitro aggregation in mutants A4V, G93A, 
G85R [63]. 

We found the dimer binding free energy of apo WT protein 
to be about 62 kJ/mol, which is in quite good agreement 
with experimental estimates derived from equilibrium de- 
naturation data on C6A/C111A double mutants by extrap- 
olating the destabilizing effects of the C6A/C111S double 
mutation (AG 2 m^m 2 ~ 61 kJ/mol [52], see Table 2). Fits 
of equilibrium denaturation data of the pseudo-WT mu- 
tant protein C6A/C111S to three-state models give some- 
what smaller numbers for dimer stability for the double 
mutant of about 51-52 kJ/mol [48,61,52]. Desideri and col- 
leagues [64] report a dimer stability for apo WT protein of 
about 52 kJ/mol, however these studies also predicted mod- 
erately smaller dimer affinity for the holo protein than the 
apo protein (about 49 kJ/mol), in contrast to our results. 

The ALS-associated mutant S134N has been observed exper- 
imentally to have only moderate effects on thermodynamic 
stability [12], indicating mechanisms aside from global un- 
folding play a role in disease propagation. Indeed, we observe 
significant effects on both metal affinity and dimer stability 



for this variant: the dimer stability is reduced from the WT 
by 45 kJ/mol, and the Zn binding free energy is reduced 
from the WT by 12 kJ/mol. Matthews and co-workers also 
see a Zn binding free energy for an obligate monomeric form 
of this mutant (C6A/C111S/F50E/G51E/S134N) that is re- 
duced from the C6A/C11S/F50E/F51E form by 7 kJ/mol, 
based on the ratio of their values [12]. On the other 
hand, some of their obligate monomeric mutants have higher 
free energy for Zn binding (e.g. their AAG values for 
A4V, G93A, L38V, and S134N are -2.9, -1.7, +2.9, and 
+6.7 kJ/mol respectively). This contrasts with the consis- 
tent decrease in Zn affinity for ALS mutants that we have 
seen here, and may be due to the mutations present in 
the their pseudo-WT background. That said, our values of 
AAG for Zn-binding free energy correlate very well with 
theirs (see Table SI): r = 0.99, P r = 0.01. Consistent with 
Crow et. al. [54], we do see reduced metal affinity for all 
fALS mutants considered, however our values for Zn bind- 
ing free energy did not correlate with their values (of A4V, 
L38V, and I113T), possibly due to monomerization during 
the course of the affinity measurements [12]. Our values 
for Cu binding free energy did correlate with their val- 
ues (r = 0.78), however the number of mutants (3) is not 
enough to be statistically significant. 

The binding free energy of Zn to WT SOD1 monomer 
as measured experimentally by Dokholyan and co-workers 
is about 40 kJ/mol [46], which is in excellent agreement 
with our value of 40 kJ/mol. However, their value was 
obtained at pH 3.5, which would reduce the value from 
that at pH 7. Likewise, calorimetry measurements at pH 
5.5 by Valentine and colleagues give comparable numbers 
of at least 45 kJ/mol for binding the first Zn ion to apo 
SOD1 homodimer [53]. In [46], these authors also observed 
that metal depletion, facilitated by dialysis against metal- 
free buffer, promoted aggregation at 30/iM and pH 3.6. We 
cannot address whether metal depletion occurs before or 
after monomerization along the aggregation pathway, and 
this does not appear to be a settled issue. As well, the 
total free energy change from holo dimer to a monomerized, 
(Cu,E) form is independent of the sequence in which these 
changes occur, so in principle both parallel pathways may 
be undertaken and would give the same net result, provided 
that that part of the misfolding/aggregation reaction is 
under thermodynamic control. 

Our calculations for dimer stability have 14 fALS mutants 
in common with the calculations of Khare et al [65], how- 
ever we see no correlation between our apo dimer stability 
values for these mutants (Table SI) and their values for 
either the free energy of dissociation of the dimer, or the 
overall thermodynamic stability of the dimer (r = 0.001 and 
r = 0.08). These authors mention that one possible source 
of discrepancy between their values and experimentally- 
determined stabilities is that their simulations relied on the 
assumption that there is no major structural change in the 
protein upon mutation. Their values also rely on calcula- 
tions of absolute free energies [66,67], for which it is often 
challenging to obtain quantitative accuracy, in contrast to 
relative changes in the free energy. We have sought to be 
careful here in calculating relative changes in free energy, in- 
cluding relaxation after mutation or PTM by thermal equi- 
libration, and also re-equilibration in the final state after 
dimerization or metal extraction. Relaxation free energies 
in demetallation generally vary mutant to mutant, and were 
more substantial for re-equilibration of disulfide-reduced or 
Zn-depleted SOD1 (Figures S3,S4). Relaxation free ener- 
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gies for monomerization were larger than for demetallation, 
and were most significant for apo(SH) SOD1, followed by 
apo(SS) and Cu,E SOD1 (Figure S5). These variants had 
significantly increased entropy in the monomeric state from 
the dimeric state because of the loss of dimer interface con- 
straints on the zinc-binding loop. 

Meiering and co-workers have found that for fALS mutants 
E100G, G93A, G85R inserted into the pseudo-WT back- 
ground C6A/C111S, dimer stability was largely unaffected 
by mutation [47]. In contrast we found significant weaken- 
ing of dimer stability due to the mutant G85R, albeit in the 
pure WT background of holo protein. Similar investigations 
from the Meiering group have shown that the same muta- 
tions in the apo dimer did weaken dimer stability, however 
the decreases in dimer stability that we have calculated 
are substantially larger than their values [48]. These dif- 
ference may be reconciled by computational studies of the 
binding free energy for fALS mutants with explicit pseudo- 
WT backgrounds (e.g. triple mutants here); mutation of 
the cysteines may modulate the dimer stability, and alter 
free energetic differences between mutants and WT. Exper- 
imental evidence suggests however that these effects may 
be modest [52,68], however the computational study is an 
interesting topic for future work. 

Several groups, including those of Oliveberg [52], O'Halloran 
and Bertini [8], and Hart and Valentine [69], have found that 
for pseudo WT SOD, either Zn binding or formation of the 
C57-C146 disulfide bond was sufficient to induce or rescue 
dimer formation, i.e. both the absence of metals (apo) and 
reduction or ablation of the disulfide bond are necessary 
to induce monomer formation, even at low concentrations 
(~ 10/iM SOD1). Sedimentation velocity analysis [69] indi- 
cates that both holo C57S and apo-SODl both sediment as 
dimers, with apo-SODl marginally more monomer-like. Con- 
sistent with these observations, we found that both metal 
depletion and disulfide reduction had significant destabi- 
lizing effects on dimer stability. The holo(SH) dimer had 
lower binding free energy than that for either the (Cu,E) 
or (E,Zn) variants, but the apo dimer had lower dimer 
binding free energy than the holo(SH) dimer; this result 
that binding of both metals (holo vs. apo) plays a larger 
role than disulfide formation in dimer affinity is consistent 
with simulations from Dokholyan's group, based on inter- 
domain native contact data [70]. Demetallation destabilized 
the disulfide-reduced dimer by about 7 kJ/mol. However, 
we found that disulfide reduction only further destabilized 
the apo dimer by about 1 kJ/mol, in apparent contrast to 
the observation that apo(SS) SOD1 elutes as a dimer but 
apo(SH) SOD1 elutes as a monomer. One possible resolu- 
tion here is the role that increased entropy might play in 
the reduced protein in rebinding a partially disrupted in- 
terface. Further work will be required on the computational 
side to reconcile these observations. 

The mutants H46R, H46R/H48Q, and D125H were atypi- 
cal in that the Zn affinity was larger than the Cu affinity 
(though only marginally for D125H). The first two of these 
variants have mutations in Cu-coordinating residues, so the 
reversal in affinity was not surprising. For D125H the rea- 
sons are more subtle however. The center of mass of the 
aspartic acid side chain is 12A from the Cu and 11 A from 
the Zn respectively, so the coupling to binding free energy 
involves bridging interactions through at least 2 side chains. 
This system provides another example of non-local commu- 
nication in the protein through stress-strain networks. 



Cu 1+ affinity is quite modest for mutants, indicating ready 
release of monovalent Cu in the reduced phase of the cat- 
alytic cycle. Released or exposed Cu is capable of react- 
ing nonspecifically with a variety of substrates to produce 
reactive oxygen and nitrogen species that are toxic to in- 
tracellular structures, including microtubules, metabolic en- 
zymes, and signalling proteins [71]. Reactive species may 
oxidatively modify side chains to inactivate SOD1 [72] and 
induce its misfolding [27]; positive feedback loops have been 
identified between protein misfolding and excitotoxicity [73]. 

By simulating AFM-like probes, we investigated the me- 
chanical stability of holo, holo(SH), apo, and apo(SH) vari- 
ants of WT protein. Here we found that the apo variant was 
less mechanically stable than holo(SH) (statistical signifi- 
cance p ~ 3 x 10 — 8 ), indicating that metal depletion is more 
mechanically destabilizing than disulfide reduction. This 
is consistent with differential scanning calorimetry (DSC) 
measurements of the melting temperatures of SOD1 vari- 
ants, where (E,Zn)(SH) SOD1 was observed to have larger 
T m than the apo variant by about 9°C [59]. Similarly, 
Cu,E(SS) is more mechanically stable than holo(SH) SOD1 
(Figure S6), and we anticipate that DSC and other mea- 
surements of this system would recapitulate the results of 
apo(SS) and holo(SH) (Table 2). Intriguingly, we found that 
apo(SH) monomer was more mechanically stable than apo 
monomer (p « 8 x 10 -8 ). This observation was supported 
by time-resolved simulations of the relaxation kinetics of 
the internal potential energy. These showed that the protein 
internal potential energy increased upon disulfide reduction 
for the holo protein, indicating decreased stability through 
loss of stabilizing interactions. However, the protein internal 
potential energy decreased upon disulfide reduction for the 
apo protein, indicating increased stability and relaxation of 
frustrating interactions. The existence of frustration in the 
apo state is supported by analysis of the mechanical stabil- 
ity of C-terminal truncation mutants, native basin dynamic 
fluctuation analysis, and frustration/potential energy anal- 
ysis [5]. The potential energy increase after SS-reduction in 
the holo protein indicates that entropy increase plays a role 
in the lowering of the free energy after the removal of con- 
straints. The potential energy decrease after SS-reduction 
in the apo protein indicates that the release of native stress 
upon removal of constraints is sufficiently large to be ob- 
served amidst the concomitant effects arising from entropy 
increase. 

These results, at first sight, contrast with DSC measure- 
ments indicating that melting temperatures are larger for 
WT apo SOD1 than for apo(SH) SOD1 by about 7°C [59], 
i.e. the disulfide bond thermodynamically stabilizes the apo 
monomer. One important difference between the simulation 
probes and experimental measurements is that our mechan- 
ical probes measure stiffness in the native basin, but not the 
relative difference in the free energies between the folded 
and unfolded states. While these generally correlate for one 
variant at different temperatures, they need not correlate 
for different protein variants at the same temperature, in 
particular if one variant has larger entropy in the unfolded 
state, as is the case here because of the presence or absence 
of a disulfide constraint. It is possible to observe a stiffer 
native basin for a less thermo-stable protein, as indicated 
schematically in inset (c) of Figure 6. 

Zinc removal dominates the effects of metal depletion on 
mechanical stability. In fact depletion of both metals results 
in lower mechanical stability than (Cu,E) SOD1 only for 
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a few (~ 8) of the least stable residues (i.e. the statistical 
significance that the cumulative distribution of apo is weaker 
than (Cu,E) is about p « 9 x 1CT 3 ). 

Insertion of metals into apo(SH) SOD1 indeed makes some 
regions very mechanically stable, increasing the stability of 
the most stable regions of the protein. However the general 
trend over most of the residues is to mechanically destabi- 
lize the protein due to added internal stresses. Thus, metal 
depletion modestly stabilizes the disulfide-reduced protein 
mechanically (p « 8 x 10 — 3 ). Together with the observation 
that disulfide reduction mechanically stabilizes apo SOD1, 
this evidence points towards cooperativity between metal 
binding and disulfide formation (Table 2). As mentioned 
above, mechanical stability need not be equivalent with ther- 
modynamic stability: thermodynamic stability accounts for 
the free energy of the unfolded state, while our mechanical 
stability measurements do not. 

For the Zinc-depleted form, shifting of the Cu to the Zn po- 
sition resulted in a more robust mechanical stability profile 
(p«lx 10~ 7 ). Supporting this trend in stability, WHAM 
simulations of the transfer of Cu from the putative Cu bind- 
ing site to the Zn binding site for the (Cu,E) form of the 
protein lowered the free energy by « —6 kJ/mol, indicating 
that the Zn binding position is more favorable for the Cu 
when only Cu is bound. This provides further evidence that 
binding of Cu is under kinetic control, modulated essentially 
by the initial Zn-coordinated folding of the protein and its 
subsequent binding to the co-expressed Cu chaperone CCS, 
which selectively loads Cu to its putative site [74]. Com- 
puting the transfer free energy for Zn from the Zn-binding 
site to the Cu-binding site gives « +6 kJ/mol, supporting 
the model that Zn-binding is thermodynamically controlled. 
Kinetics experiments by Oliveberg and colleagues indicate 
Zn 2+ binds to the putative Cu site during folding, but may 
be transferred to the higher affinity Zn-binding site late in 
the folding process [75]. We note that the above transfer 
free energies are estimates, and do not take into account 
the differential affinity due to the 3d 10 vs. 3d 9 electronic 
valence of Zn 2 + and Cu 2+ respectively. 

The decrease in free energy corresponding to the shift of the 
Cu is accompanied by increased mechanical rigidity in the 
Zn binding loop. However, some residues, such as residue 
24, 75 and 120, are significantly weakened by the shift. Con- 
sistent with significantly weakened mechanical rigidity, the 
root mean squared fluctuation (RMSF) values are increased 
for those residues, corresponding to an increase in local en- 
tropy. Similar effects were seen here for processes such as 
metal binding, and disulfide formation- energetically favor- 
able processes that result in local enhancement of struc- 
ture, but which may increase entropy non-locally. Non-local 
entropy generation has also been seen in pulling simula- 
tions, where local mechanical stress induced by pulling on 
a particular residue resulted in non-local unwinding of the 
short stretches of a helix in the Zn-binding and electrostatic 
loops [5]. 

These effects perhaps foreshadow one of the essential fea- 
tures of disease propagation in template-directed misfolding, 
namely that of non-local entropy transduction across the 
protein. It is feasible that the induced misfolding of a pro- 
tein due to the interaction with misfolded template would 
correspond to low enthalpy and low entropy in the vicinity 
of the binding interface, and high enthalpy and high en- 
tropy away from the binding interface. Such a phenomenon 



would facilitate the structural plasticity required to ren- 
der the protein amenable to adopt altered conformations in 
the presence of misfolded template. Once altered conforma- 
tions are adopted, the misfolded protein becomes part of 
the infectious species, capable of further seeding misfolding 
through the reservoir of healthy protein. 



4 Methods 



Jj. . 1 Steered molecular dynamics simulations 



Steered molecular dynamics (SMD) with constant-velocity 
moving restraints involving two tethering points was used 
to simulate the action of a moving AFM cantilever on a 
protein. Three computational assays were investigated: 

1. A study of SOD1 monomer, where one tether was placed 
at the position of the alpha carbon closest to the center of 
mass of the protein (generally C a (46)), and another tether 
was placed on the alpha carbon of a particular amino acid 
to be pulled on. 

2. A study of monomerization from the dimer, with a tether 
on the alpha carbon atom closest to center of mass of each 
monomer. 

3. A study of metal removal, with one tether on the C Q 
atom of either residue Phe45 (for Cu removal) or that of 
Asp83 (for Zn removal), and the other tether on either the 
Cu or Zn metal ion. These residues are chosen because 
they determine the approximate direction of highest solvent 
exposure of the metal. 

In assay (1), to implement a mechanical stability scan across 
the surface of the protein, every 5th amino acid along 
with the first and last residues was selected for a pulling 
simulation. This coarse-grains the mechanical profile, so that 
the work values obtained are a sampling of the work values 
for all residues. 

As described in [5], SMD simulations [76-79] were preformed 
in GROMACS using an all-atom representation of the pro- 
tein, with OPLS-AA/L parameters [80,81], and a general- 
ized Born surface area (GBSA) implicit solvent model with 
protein dielectric constant 4, solvent dielectric constant 80, 
and surface tension coefficient for solvation free energies of 
0.005 kcal/mol/(A 2 ) [82]. Born rad ii are calculated from the 
Onufriev-Bashford-Case (OBC) algorithm [83], and bond 
lengths containing a hydrogen atom are constrained by the 
LINCS algorithm [84]. The long-distance cut-off used for 
non-bonded interactions was 14 A for both Electrostatic 
and van der Waals (VDW) interactions. 

After energy minimization by steepest descent to remove 
potential steric clashes, simulations were implemented using 
an integration time-step of 2 fs, in a cubic simulation box 
with periodic boundary conditions and size such that all 
protein atoms were initially at least 20 A from any cubic 
face. An equilibrium simulation was always run for 20 ns to 
equilibrate the protein before any pulling simulations were 
performed and data was collected; coordinates were saved 
every 100 ps. 

Because only moderate deformations were required to con- 
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struct mechanical force-extension curves, a very slow pulling 
speed (for simulations) of 2.5 x 10 — 3 m/s was admissible 
when work values were computed. This speed is 400 — 8000 x 
slower than the speeds generally used in simulations [85- 
87], and is comparable to pulling speeds used in AFM ex- 
periments [88]. The spring constant of the simulated AFM 
cantilever was taken to be 5 kJ/mol/A 2 , which is compa- 
rable to the protein effective spring constant. 

For assay (1), each pulling simulation ran until the change 
in distance between the two tethering points was 5 A. For 
assays 2 and 3, the simulations were continued until the 
force fell to a value of zero indicating the separation of the 
dimer or demetallation of the monomer, respectively. 



Jf.,2 Tethering residues 



The Cot atom of residue 46 was closest to the center of 
mass for WT SOD1 and all variants. For mechanical scans, 
residues 40, 45, 50, 54, and 55 were either too close to 
the central C a , or were along the same beta strand and 
gave anomalous force-extension profiles with artificially large 
forces probing covalent bonding topology moreso than non- 
covalent stabilizing interactions. In these cases, the tethering 
residue was moved to residue 76. This always resolved the 
problem of large forces. In assay (2) for dimer pulling 
simulations, the tethering point (on each of the monomers) 
was C a (46). 

In assay (3) for Cu and Zn pulling simulations, the tethering 
point was C a (45) for Cu extraction and C a (83) for Zn. 
The residue selection was chosen by finding residues within 
5A of the metal, and selecting the residue that was the 
most buried (least solvent accessible surface area (SASA)) 
as the tethering residue. The resulting pathways are shown 
in the insets to Figure 3 (a). For the analysis of the (Cu,E) 
variant PDB 2R27, the Cu ion of chain A was shifted by 
1.3 A to a location associated with reduced (Cu 1+ ) that 
is closer to the putative Zn position. For this protein the 
tethering residue was taken to be C a (83). 



Proteins considered 



Proteins whose mechanical profiles were calculated are given 
in Table 1 below. For (Cu,E) SOD1, 2R27 was used to 
generate an initial structure with the Cu shifted to the Zn 
position (as it exists in the PDB structure), and a (Cu,E) 
variant with the Cu in its putative position, by moving the 
Cu and equilibrating the structure. 

Several mutants were also considered for the calculation 
of metal affinity and dimer affinity. These included A4V, 
D124V, D125H, D76Y, D90A, G127X, G37R, G41D, G41S, 
G85R, G93A, G93C, H43R, H46R, H46R/H48Q, H80R, 
I113T, L144F, L38V, S134N, T54R and W32S. The corre- 
sponding PDB entries used to generate structures for these 
mutants for simulations are given in Table 1. 



^..5 Remodeling disordered/missing regions in PDB 
structures 



(Cu,E) SOD1 (2R27) has two unstructured segments in the 
crystal structure, residues 68-78 and residues 132-139, how- 
ever a full length construct must be used as an initial condi- 
tion in simulations. The missing segments were added to the 
protein by the following method: First, segments of protein 
corresponding to the missing sequences were taken from the 
native holo SOD1 structure, and Replica Exchange Molec- 
ular Dynamics simulations of the free peptides in explicit 
solvent [90] were performed using the NAMD simulation 
package [91], to obtain an ensemble of disordered conforma- 
tions for these two peptide fragments. Clustering was then 
performed on the simulated ensemble for each peptide [92], 
and a conformation was chosen from the largest cluster, such 
that its end-to-end distance best matched the corresponding 
distance between the structured residues present in 2R27 
that bracketed the missing sequence. After adding the miss- 
ing loop segments, we "back- mutated" four residues, A6C, 
S80H, S83D, and S111C, to recover the WT sequence, using 
the PyMOL software package [93]. The structure was then 
energy minimized, and equilibrated for 20 ns, as described 
in section 4.1. A full length construct was thus assembled 
to be used as an initial condition in pulling simulations. 



J[.3 Residues used for mechanical scans of protein sta- 
bility 



Every 5th residue along the protein sequence, including the 
first and last residues, was taken as a tethering residue. 
As well, the central positions of the predicted epitopes and 
anti-epitopes were added to the data set [89], as described 
in reference [5]. If the center of the epitope happened to 
coincide with a multiple of 5 already included in the scan 
(e.g. residue 10), then the next residue (residue 11 in this 
case) was also taken. Thus a set of 48 residues was used in 
the mechanical scans: (1, 5, 10, 11, 15, 17, 20, 24, 25, 30, 
31, 35, 38, 40, 45, 46, 50, 54, 55, 60, 65, 70, 73, 75, 80, 85, 
90, 91, 95, 100, 101, 105, 108, 110, 115, 120, 121, 125, 130, 
135, 136, 140, 141, 145, 146, 150, 151, 153). 



Jj. . 6 Modeling proteins with no PDB structure 



Some proteins studied here had no PDB structural coordi- 
nates available; see Table 1. For these proteins, a structure 
was built by modifying known PDB structures of similar pro- 
teins. For example, disulfide-reduced WT SOD1 was created 
from (E,Zn)(SH) SOD1 (2AF2) by adding Cu to that struc- 
ture, and back-mutating A6C and S111C using the PyMOL 
software package [93]. Apo WT SOD1 was created from the 
apo SOD1 structure PDB 1RK7 by back-mutating Q133E, 
E50F and E51G. Apo(SH) SOD1 was created from the apo 
WT SOD1 structure 1RK7 by reducing the disulphide link- 
age and mutating E50F, E51G, and Q133E. The mutant 
H46R was generated from the H46R/H48Q double-mutant 
of SOD1 by back-mutating Q48H. Similarly, the mutants 
W32S, G41S, G41D, G93C, L144F, D90A and D76Y were 
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prepared from the hole- SOD1 by mutating at the respec- 
tive positions. The G127X mutant was created from holo 
SOD1 by deleting the last 20 residues, and then mutating 
the last 6 residues of the remaining sequence (KGGNEE) 
to correspond to the frame-shifted non-native C-terminal 
peptide sequence (GGQRWK) prior to the termination se- 
quence. All modified structures were energy minimized and 
equilibrated for 20 ns by running an equilibrium simulation, 
before performing any pulling assays. 



Jj.,7 Umbrella Sampling and Weighted Histogram Analy- 
sis Method (WHAM) 



Jj.,7.1 WHAM procedure for mechanical profiles 

Initial configurations for use in WHAM were obtained from 
the pulling simulations described in Methods section 4.1: 25 
initial conditions between and 5A were simulated for 10ns 
each, in an umbrella potential with stiffness 500kJ/mol/nm 2 
to constrain the simulations near their corresponding sep- 
aration distances. These 25 simulations are then used to 
reconstruct the free energy profile along the distance co- 
ordinate using the weighted histogram analysis method 
(WHAM) [94], and the free energy difference between 
and 5A is then obtained. A convergence check assured the 
result was unchanged if 30, 35 or 40 windows were used. 



Jf.,7.2 WHAM procedure for metal expulsion and dimer 
separation 

To obtain free energies of metal binding, metals were in- 
serted into the putative binding locations for all SOD1 
variants, where they were found to be at least metastable. 
Dimer stabilities were calculated under the same metallation 
conditions as the PDB structure for all mutants. Mutants 
with no structure in the PDB were fully metallated. 

For the metal expulsion free energy calculations in Sec- 
tion 2.2, a most-direct straight-line path was first deter- 
mined for the pulling direction as described in Methods 
section 4.1, which determined the tethering residue. The 
metal was pulled away from the protein using a spring con- 
stant of 500 kJ/mol/nm 2 and a pull rate of 0.01 nm/ps. For 
purposes of the free energy calculation, the final extension 
from the equilibrium distance between metal and tethering 
residue was taken to be approximately 30 A. From these 
trajectories, snapshots were taken to generate the start- 
ing configurations for the umbrella sampling windows [95- 
97]. An asymmetric distribution of sampling windows was 
used, such that the window spacing was lA between and 
20 A separation, and 2A beyond 20A of separation. Such 
spacing allowed for increasing detail at smaller separation 
distances, and resulted in a total of 25 windows. In each 
window, 10 ns of MD was performed for a total simula- 
tion time of 250 ns utilized for umbrella sampling. Analysis 
of the results was performed using the weighted histogram 
analysis method (WHAM) [94]. 

To remove conformational distortion effects on the free en- 
ergy, position restraints were applied in the following way. 
For a given C a atom in the protein, all other C a atoms 



within 5 A were constrained to have a roughly constant C a - 
Ca. distance, the same as in the equilibrated structure, us- 
ing spring constants of 392 x 10 3 kJ/mol/nm 2 . The spring 
constant was varied to interpolate between regimes where 
protein deformation is dominant, and a regime where the 
constraining C a -C a network can potentially influence the 
expulsion free energy. Using C a constraints allows the pro- 
tein to retain its structure under force, while still allowing 
the side chains and partially the backbone to fluctuate in 
response to the external perturbation. After metal expul- 
sion, Co. constraints were relaxed and the relaxation free 
energy was calculated as described below. 

A similar procedure was employed for dimer separation. 
Here the tethering points were taken to be the C a atom 
of residue 46 in each of the monomers. The monomers 
constituting the dimer were pulled away from each other, 
until the final distance was approximately 30 A. The spring 
constant, pull rate, and window separation were the same 
as in the metal expulsion analysis. Within each monomer, 
C a constraints were also applied to fix the relative positions 
of the C a atoms near their equilibrium positions, and thus 
minimize conformational distortion during monomerization. 

The free energy for metal expulsion or dimer separation 
was corrected to account for protein relaxation in the final 
metal depleted or monomerized state: AFtot = AF CO nstr + 
AF re i ax . That is, after either metal expulsion or monomer- 
ization, Col constraints are gradually reduced from 392 x 10 3 
kJ/mol/nm 2 to zero using 30 windows and 10ns of relax- 
ation time in each window, and free energy changes for this 
process are again obtained using WHAM. Convergence of 
the free energy values was tested by both varying the num- 
ber of windows (to 30, 35 and 45) and varying the length 
of the equilibrium simulation in each window (to 15, 20, 25 
and 30 ns). In all cases the free energies were seen to have 
converged using the original protocol. 

To further confirm that equilibrium has been reached, ther- 
modynamic cycles were constructed for either re-insertion of 
the metal, or re-dimerization of the monomers. The metal or 
dimer binding free energies subject to C a constraints were 
calculated, as well as the subsequent relaxation free ener- 
gies. Schematics of the thermodynamic cycles are plotted 
in Figure 4, and values of each of the free energy changes 
for WT protein and several mutants are given in Table S2 
of the Supplementary Content. 



WHAM procedure for shifting Cu to the Zn- 
binding position 

In the simulated conformations of (Cu,E) SOD1 that have 
been equilibrated from the crystal structure 2R27 (in which 
Cu is bound close to the putative Zn position), Cu remains 
close to the Zn position and does not shift back to the Cu 
binding site. The distance between the position of the Cu ion 
in the equilibrated structures from simulations of holo WT 
SOD1 and the above (Cu,E) SOD1 is 7.92 A. We calculated 
the free energy change to move Cu from its putative binding 
position to the Zn-position with the following method. We 
first pulled the Cu from its original position to the (Zn- 
binding) equilibrated position by applying a spring constant 
of 500 kJ/mol/nm 2 and a pulling speed of 10 m/s, while 
enforcing harmonic constraints on all of the C a atoms as 
in Methods Section 4.7.2. Tethers were placed on the Cu 
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ion and the C a atom of residue 57, because it is almost 
on the line joining the initial and final positions of the Cu. 
Conformations were collected after every 0.5 A, and for each 
of the conformations a 10 ns MD simulation was performed. 
These trajectories were used for umbrella sampling, and 
analysis of the results were performed using WHAM. 



Jj. . 8 Free energy differences from Jarzynski equalities, and 
finite sample size corrections 



The free energy cost AF to extend a residue to 5A is found 
from Jarzynski 's equality 



Jj. . 9 Statistical Analysis 



To compare whether the difference between two distinct cu- 
mulative distributions is statistically signicant, the method 
described in reference [5] is used. A summary of the method 
follows. The mechanical profile is a collection of 48 work 
values for a given SOD1 variant, where each work value has 
an error of a = 2.7kJ/mol. We wish to answer the following 
question: given one work profile and its resultant cumulative 
distribution, along with the error bars associated with the 
individual work values, what is the probability of obtaining 
a cumulative distributions at least as extreme as another 
given one by chance, i.e. as arising from the original work 
profile? 



AF = -kT\n(e- w ' kT ) , 



where W is the work from a non-equilibrium measurement, 
kT is Boltzmann's constant times the temperature in Kelvin, 
and the average (• • •) is the ensemble average over replicated 
pulling assays. For a finite sample size, the free energy in 
Eq. (1) becomes 



AF, 
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-W/kT 



,(2) 



which must be corrected because of systematic bias for finite 
sample size. Following Gore, Ritort, and Bustamante [42], 
an estimate for the free energy change AF that accounts 
for finite sample-size corrections is given by the following 
set of equations: 



AF = AF T - 



W dl 



is2 



w 



w d 



, W di8 



N ^( W dis) 

W dis = (W) N - AFj 

a (W) = In (2W/kT) / In [C (e 2 w l kT - l) ] 



(1) To find the above probability and thus determine the sta- 

tistical significance of a given cumulative distribution, we 
employ a "Monte Carlo" procedure of generating cumula- 
tive distributions from a given "baseline" work profile. For 
example, to test the hypothesis that shifting the Cu from its 
putative position to the Zn-binding position mechanically 
stabilizes the (Cu,E) protein, we wish to find the probability 
of obtaining a cumulative distribution at least as stabilized 
as the (Cu-shift,E) cumulative distribution from the original 
(Cu,E) cumulative distribution. The (Cu-shift,E) cumulative 
distribution has 31 work values that are more stable than 
those in the (Cu,E) cumulative distribution; the difference 
in work values are plotted in Figure S2c, largest to smallest. 
We seek cumulative distributions that have rank ordered 
deviations at least as large as those of the (Cu-shift,E) 
cumulative distribution. Work profiles are constructed by 
adding gaussian noise of mean zero and standard devia- 
tion a = 2.7 kJ/mol to a the baseline work profile (here 
that of (Cu,E) SOD1). A cumulative distribution is then 
constructed for each generated profile by sorting the values 
lowest to highest. After constructing N of such cumulative 
distributions, we ask whether one has found a cumulative 
distribution with, e.g. rank-ordered positive deviations at 
least as large as those in the (Cu-shift,E) cumulative dis- 
tribution. The value of N where the expected number of 
trajectories is unity determines the statistical significance 
(3) of the difference between the two profiles: p = 1/N. Plots of 

the corresponding distributions for (Cu,E) SOD1 and (Cu- 
shift,E) SOD1 are given in Figure S2 of the Supplementary 
Content. 



In the above equations, {W) N is the average work for the 
N pulling assays, and Wdis is the average dissipated work- 
work that did not go into producing the free energy change. 
The function a(W) is evaluated at two values of dissipated 
work, Wdis, an d a corrected estimate for the dissipated 
work, Wdis2, which accounts for the fact that Wdis is 
itself an underestimate because of the finite sample size 
bias. The constant appearing in a(W) defines a crossover 
between small N and large N regimes and is not rigorously 
determined. We take C = 15 following Gore et al [42]. 

Table S6 in the Supplementary Content gives numerical 
values for the various quantities in Equation (3), in the 
assay where residues 10 and 17 are pulled to 5 A. The 
mean dissipated work W ' dis2 m our pulling assay is only 
about a kJ/mol, indicating that pulling is near equilibrium, 
and finite sample-size corrections are not large: about 1-3 
kJ/mol. 
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Fig. SI. The effect of metal depletion on the mechanical stability of SOD1. (Inset panel a) Mechanical profiles of holo SOD1, (E,Zn) 
SOD1, (Cu,E) SOD1, and Zn-depleted SOD1 with the Cu ion shifted from its putative position to that of the Zn (Cu-shift,E), which 
is near the location of the Cu observed in the crystal structure of PDB 2R27. Profiles are color-coded according to the legend in the 
upper left. The effect of metal depletion modulates the entire mechanical stability profile. Work values for individual residues are given 
in Table S5 of the Supplementary Content. (Inset b) Probability distributions of mechanical stability for holo and (Cu,E) variants. This 
representation shows that the distributions are indeed different, however there is significant overlap, which along with the necessity of 
binning the data to construct the histogram, makes the differences difficult to quantify. (Main Panel) Cumulative histograms of the 
above SOD1 variants. The cumulative histogram requires no binning, and simply rank orders and accumulates the work values. Using 
this representation, discrepancies between work profiles most clearly emerge. The mechanical stability of the SOD1 variants can be rank 
ordered, strongest to weakest, as holo, (E,Zn), (Cu-shift,E), and (Cu,E). The distribution of apo SOD1 is broadened compared to (Cu,E) 
SOD1, but not overall weakened. Also shown in the main panel is a correlation matrix resulting from a pairwise comparison of the work 
values for all pairings of SOD1 variants. The profiles do not correlate, again supporting the notion that metal depletion modulates the 
entire mechanical stability profile. 
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Fig. S2. (a) Statistical significance test between (Cu,E) and (Cu-shift,E) cumulative distributions. The black curve denotes the 
(Cu,E) profile, and the red curve denotes the (Cu-shift,E) profile. The collection of green lines corresponds to 100 randomly-generated 
cumulative distributions in the following way: starting from the (Cu,E) work values, gaussian noise with zero mean and standard deviation 
a = 2.7 kJ/mol is added to each work value. The resulting collection of work values are then rank-ordered to generate a cumulative 
distribution, and this process is then repeated TV times. The collection of blue lines corresponds to TV = 1000. For large enough N, one 
will begin to find outlying cumulative distributions which deviate as much or more than a given trial distribution, here (Cu-shift,E). An 
example of such an outlier is shown in panel (b). In this case, the number of positive deviations (differences in work value) between the red 
and black curves in panel (a) are recorded (31), along with their work values, which are themselves rank ordered largest to smallest and 
plotted in figure (c). An outlier at least as extreme as (Cu-shift,E) must have a set of positive work deviations such that the corresponding 
curve lies above the red curve in panel (c). In this case the randomly generated outlier has 36 positive deviations in work, 31 of which 
are larger than those between (Cu-shift,E) and (Cu,E) variants. The number of randomly generated trajectories is increased, until the 
mean number of successful trajectories, corresponding to the criterion in panel (c), exceeds unity. A plot of the mean number of successful 
trajectories as a function of TV is given in panel (d). The statistical significance p of a given trajectory is then given as p = 1/TVi, where Ni 
is the number of randomly generated trajectories that give an expectation value of unity for the number of successfully-generated outliers. 
In this case, the statistical significance of (Cu-shift,E) is about 1/(7 x 10 6 ), or about 1.4 x 10 -7 . 




Fig. S3. Free energy change for Cu expulsion, for WT SOD1 along with several mutants and post-translationally modified variants. The 
heights of the blue bars indicate the corresponding free energies of metal binding before conformational relaxation of the protein in the 
unbound state. Red bar heights indicate the corresponding free energies of binding after relaxation of the protein in the unbound state. 
A decrease in free energy value is observed for all the proteins. 




Fig. S4. Free energy change for Zn expulsion, for WT SOD1 along with several mutants and post-translationally modified variants. The 
heights of the blue bars indicate the corresponding free energies of metal binding before conformational relaxation of the protein in the 
unbound state. Red bar heights indicate the corresponding free energies of binding after relaxation of the protein in the unbound state. 
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Fig. S5. Free energy change for dimer separation, for WT holo SOD1 along with several mutants and post-translationally modified 
variants. All mutants plotted here are taken in the holo state. The heights of the blue bars indicate the corresponding free energy to 
monomerize the dimer before conformational relaxation of both proteins in the unbound state. Red bar heights indicate the corresponding 
free energies of binding after conformational relaxation of both protein in the monomerized state. 
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Fig. S6. Loss of Zn is at least as destabilizing to SOD1 mechanical stability as disulfide reduction. Panel (a) gives the cumulative 
distribution of the work profiles for holo(SH) and Cu,E(SS) SOD1 variants, which shows that Cu,E(SS) SOD1 has marginally weaker 
mechanical stability than holo(SH) SOD1 (p = 1.4e-7). This trend is also observed in the cumulative distribution of the equilibrium 
fluctuations of the holo(SH) and Cu,E(SS) variants (Panel (c)), where Cu,E(SS) shows marginally larger fluctuations. Panel (b) compares 
cumulative distributions of work profiles for holo(SH) and apo(SS) SOD1, showing apo(SS) SOD1 is marginally weaker mechanically 
(p = 3e-8). The cumulative distribution of equilibrium fluctuations of these two variants (panel d) also shows that apo(SS) SOD1 has 
more significant dynamics in the native state. 
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Fig. S7. Zn-binding and electrostatic loops (ZBL and ESL) show decreased mechanical rigidity and increased dynamic mobility upon 
metal release. The SOD1 sequence is divided into 5 different regions: the region N-terminal to the ZBL (residues 1-48), the ZBL (residues 
49-83), the middle region of sequence between the ZBL and ESL (residues 84-120), the ESL (residues 121-142), and the C-terminal region 
after the ESL (residues 143-153). Panel (a) gives the average work values (as calculated by pulling to 5A) for those regions, for holo, (E,Zn), 
(Cu,E) and apo SOD1 variants. Panel (b) gives the difference in average work values of those 5 regions in the various metal-deficient 



forms of SOD1 from holo SOD1 (e.g. (W) h 



(W) a ). This shows that metal loss systematically destabilizes the ZBL and ESL; other 



regions of the proteins show less systematic trends. Panel (c) gives the root mean squared fluctuations (RMSF) averaged over the above 5 
regions, and panel (d) gives the difference in RMSF of those 5 regions in the various metal-deficient forms of SOD1 from holo SOD1 (e.g. 
(RMSF) apo — (RMSF) holo ). This recapitulates the mechanical stability assay: metal loss systematically increases dynamic fluctuations 
in the ZBL and ESL. Other regions show weaker trends in dynamic fluctuations. 



Supplemental Tables 



Table SI. Free energy of Cu, Zn expulsion and 
Dimer separation (kJ/mol) 
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Table S2. Free energy change for metal expulsion and dimer 
separation (kJ/mol) - thermodynamic cycle (see Figure 4, 
main text) 



Metal 


j^iconstr 




-piconstr 




WT(Cu) 


55.0 


2.1 


53.5 


0.8 


WT(Zn) 


45.1 


5.5 


41.4 


1.9 


Cu,E(Cu) 


45.6 


4.1 


43.4 


1.7 


E,Zn(Zn) 


40.2 


6.5 


36.6 


3.0 


A4V(Cu) 


46.9 


2.1 


45.5 


0.9 


A4V(Zn) 


41.9 


5.5 


38.6 


2.1 


G93A(Cu) 


51.8 


2.1 


50.5 


0.9 


G93A(Zn) 


40.6 


5.6 


37.4 


2.2 


Dimer 


^jpconstr 


^pretax 


^jpconstr 


_ AF relax 


WT 


120.1 


15.2 


111.7 


6.6 


Cu,E 


95.8 


19.2 


86.0 


9.2 


E,Zn 


110.2 


17.1 


100.2 


7.2 


A4V 


91.0 


15.3 


84.9 


8.9 


G93A 


102.1 


15.3 


93.7 


7.1 



Table S3. Correlation between Work and RMSF of SOD1 variants 



Correlation 


holo 


E,Zn 


Cu-shift,E 


Cu,E 


apo 


holo(SH) 


apo(SH) 


r 
P 


0.06 
0.69 


-0.02 
0.89 


0.06 
0.66 


0.08 
0.57 


0.01 
0.94 


0.07 
0.61 


4.7 x 10" 5 
0.99 



Table S4. Correlation table of work values for different variants of SOD1 



\. rt 
P ^\ 


holo 


E,Zn 


Cu-shift,E 


Cu,E 


apo 


holo(SH) 


apo(SH) 


holo 


1 

1 


0.25 


0.32 


0.21 


-0.37 
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-0.16 


E,Zn 


0.09 


1 


0.25 


0.15 


0.16 
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-0.07 


Cu-shift,E 
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0.09 
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0.0002 


0.02 
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'r indicates the correlation coefficient between work profiles of different SOD1 variants, and p indicates the statistical significance of the correlation coefficient. 



Table S5. Mechanical work values of WT S0D1 variants - the values reported 
in the table depict the work in kJ/mol needed to pull a particular residue up 
to 5 A. 



Res.lnd. 
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94.6 


79.4 


81.2 


11 


108.4 


98.9 
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66.5 
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87.9 
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68.1 


80.1 


103.1 
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79.4 


116.2 


64.1 


55.8 


116.0 


104.2 


85.1 
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62.5 


88.8 


54.4 


83.6 


46.7 


61.0 


70.2 
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57.3 


87.9 


84.1 


77.4 


79.9 


92.1 
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147.3 


98.9 


66.6 


76.3 


58.3 


54.8 


56.8 


145 


72.8 


84.4 


54.3 


51.7 


62.2 


71.2 


76.2 


146 


80.5 


107.0 


92.4 


59.3 


45.3 


61.0 


61.2 


150 


68.1 


55.3 


103.0 


63.3 


74.8 


99.2 


80.2 


151 


76.6 


68.0 


60.6 


54.7 


61.7 


58.0 


66.1 


153 


69.0 


52.1 


68.2 


62.1 


69.5 


60.4 


60.3 



Table S6. Finite sample size correction of Jarzynski equality (energies are in 
kJ/mol) 



Residue Index 


(W) N 


W d i S 


W d is2 


a(W dls ) 


OL{Wdis2) 




AF 


10 


56.24 


0.110 


1.33 


-0.804 


0.25 


56.13 


55.50 


17 


80.27 


0.108 


1.13 


-0.756 


0.28 


80.16 


77.55 



